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Abstract 

Infinite  periodic  structures  have  been  studied  heavily  because  of  their  efficient  filter¬ 
ing  capabilities.  They  generally  exhibit  sharp  frequency  roll-offs  at  the  frequency  band  of 
interest.  In  the  RF  region  of  the  electromagnetic  spectrum,  periodic  structures  find  appli¬ 
cations  such  as  radomes  where  they  allow  a  certain  band  of  frequencies  to  pass  through, 
andphotonic  bandgap  materials  that  block  transmission  at  selected  frequency  bands.  Most 
studies  have  been  done  with  infinitely  periodic  arrays  because  it  is  convenient  to  collapse  an 
infinite  array  into  one  representative  period  using  Floquet  Analysis.  In  reality,  infinite  arrays 
are  not  physically  realizable.  However,  truncating  an  infinite  array  introduces  an  edge  and 
invalidates  Floquet  analysis  over  the  entire  array. 

This  thesis  formulates  a  Finite  Element  Method  (FEM)  solution  of  a  semi-infinite 
periodic  array  consisting  of  infinitely  long  cylinders.  The  array  elements  sufficiently  far 
from  the  edge  are  implemented  using  the  concept  of  a  Physical  Basis  Function  (PBF).  The 
PBF  concept  is  based  on  an  a  priori  knowledge  that  the  amplitudes  of  the  currents  in  the 
periodic  elements  that  are  sufficiently  far  from  an  edge  are  constant.  Implementation  of  the 
PBF  concept  allows  the  solution  domain  of  the  FEM  to  be  bounded  by  introducing  a  periodic 
boundary  that  represents  the  truncated  portion  of  the  periodic  array. 

The  periodic  boundary  is  implemented  by  relating  the  fields  there  with  a  Floquet  phase 
factor  based  on  one  periodic  element  external  to  the  FEM  domain.  Performance  of  the  pe¬ 
riodic  boundary  at  normal  incidence  is  promising.  At  off-normal  incidence,  the  performance 
of  the  implemented  boundary  is  below  expectations.  Implementation  of  a  periodic  boundary 
by  relating  the  fields  there  with  a  Floquet  phase  factor  with  one  interior  periodic  element  is 
the  next  stage  in  the  pursuit  of  improving  off-  normal  incidence  performance. 
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APPLICATION  OF  THE  FINITE  ELEMENT  METHOD  TO  THE 
SCATTERING  OF  A  TWO-DIMENSIONAL,  SEMI-INFINITE 

PERIODIC  STRUCTURE 


/.  Introduction 

Frequency  selective  surface  (FSS)  design  and  analysis,  especially  in  the  microwave  re¬ 
gion  of  the  electromagnetic  spectrum,  has  been  studied  heavily  during  the  last  three  decades. 
Interest  in  the  exploitation  of  periodic  structures  technology  occured  because  FSSs  exhibit 
very  desirable  filtering  characteristics.  For  example,  enclosing  an  electromagnetic  scatterer 
or  radiator  within  an  FSS  provides  the  means  for  restricting  transmission  in  and  out  of  the 
enclosure  to  a  sharply  defined  band  of  frequencies.  This  effectively  “hides”  an  active  ra¬ 
diator  and  the  cavity  enclosed  within  the  FSS  from  being  detected  at  any  other  frequency 
outside  of  the  designed  transmission  band.  Shaping  the  FSS  to  conform  to  its  environment 
also  helps  minimize  specular  reflections.  Thus,  aircraft  designers  are  given  more  options  in 
developing  low  radar  cross  section  (RCS)  platforms.  These  important  FSS  characteristics 
provide  avenues  for  the  implementation  of  the  modern  version  of  an  ancient  military  strategy 
of  attacking  the  enemy  without  warning  or  very  minimal  indication  of  an  impending  attack. 
This  effective  way  of  concealing  objects  behind  a  periodic  structure  is  part  of  the  modern 
warfare  strategy  we  now  call  “stealth”  technologies. 

1.1  Background 

Although  the  science  of  periodic  structures  can  be  traced  back  to  the  end  of  the  last 
century,  accurate  calculations  were  not  achievable  until  the  introduction  of  digital  comput¬ 
ers  [20].  Faster  calculations  thus  enabled  the  deeper  study  of  these  structures  and  imple¬ 
mentation  of  new  technologies  such  as  complex  radomes,  photonic  bandgap  materials  and 
phased  array  structures,  to  name  a  few. 
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Figure  1.1  Typical  Response  of  a  Frequency  Selective  Surface  [20]. 

Frequency  selective  surfaces  behave  as  filters  of  electromagnetic  waves.  Their  transmis¬ 
sion  and  reflection  properties  vary  with  frequencies  and  incident  angles.  A  typical  response 
of  an  FSS  is  shown  in  Figure  1.1.  In  the  figure,  it  shows  that  at  a  certain  band  of  fre¬ 
quencies,  the  FSS  acts  as  either  a  bandpass  or  a  bandstop  filter,  depending  on  the  physical 
characteristics  of  its  periodic  elements. 

Periodic  structures  can  be  considered  as  one  of  two  types.  The  first  type  consists 
of  apertures  of  arbitrary  shape  “punched  out”  from  a  conducting  sheet.  The  second  type 
consists  of  conducting  elements  of  arbitrary  shape  sandwiched  inside  a  dielectric  slab.  Finite 
length  dipoles  are  the  simplest  example  for  the  second  type  [11]. 

The  usual  way  of  analyzing  periodic  structures  is  to  take  a  representative  period  and 
designate  it  as  a  reference  cell.  For  plane  wave  incidence,  Floquet’s  periodicity  condition 
stipulates  that  the  amplitude  of  the  field  at  any  point  in  the  array  is  equal  to  the  field 
in  the  reference  cell  except  for  a  phase  factor.  Essentially,  Floquet  analysis  reduces  an 
infinite  domain  problem  to  a  finite  domain  one.  Numerical  techniques  can  then  be  utilized 
to  calculate  the  structure’s  transmission 
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Traditionally,  most  “open  region”  or  unbounded  domain  problems  are  numerically 
solved  using  the  method  of  moments  (MoM).  The  moment  method  provides  a  logical  first 
step  in  the  analysis  of  the  scattering  of  periodic  structures  because  the  MoM  is  capable 
of  determining  the  scattered  fields  by  using  an  integral  equation  containing  the  free  space 
Green’s  function  that  satisfies  Sommerfeld’s  radiation  boundary  conditions.  Integral  equa¬ 
tion  techniques  in  their  purest  form,  like  the  MoM  for  instance,  are  often  restricted  to 
computing  the  scattering  from  a  limited  number  of  geometries  because  implementation  of 
the  free  space  Green’s  function  is  often  restricted  to  simple  canonical  forms.  That  makes 
the  MoM  computationally  intensive  when  calculating  complex  geometries  and  geometries 
immersed  in  non-homogeneous  dielectric  materials  without  resorting  to  hybrid  techniques. 

More  recently  in  the  computational  electromagnetics  (GEM)  field,  partial  differential 
equation  (PDE)  methods  have  come  into  greater  use  in  scattering  problems  because  of  their 
ability  to  account  for  more  complex  geometries  and  inhomogeneous  dielectric  materials  [10, 
21].  The  finite  element  method  discretizes  the  problem  domain  into  finite  cells,  also  called 
grid  or  “mesh”  elements.  These  mesh  elements  may  be  assembled  from  non-overlapping 
irregularly  sized  triangles  or  quadrilateral  elements  for  two-dimensional  space;  tetrahedra 
or  rectangular  “bricks”  for  three-  dimensional  space.  On  the  other  hand,  finite  difference 
methods  normally  discretizes  the  problem  domain  with  a  uniformly  spaced  grid,  preferably 
Cartesian  [10]. 

PDE  methods  however  require  their  problem  domain  to  be  finite.  In  order  to  character¬ 
ize  an  infinite  region  normally  associated  with  radiation  and  scattering  problems,  the  PDE 
methods  implement  a  fictitious  boundary.  Mathematically,  this  boundary  allows  incident 
energy  to  penetrate  at  the  same  time  maintaining  the  the  scattered  field  characteristics  of 
outward  propagation  to  infinity.  The  boundary  must  also  prevent  reflections  back  into  the 
solution  domain. 

In  real  world  applications,  it  is  not  possible  to  have  infinitely  long  periodic  surfaces.  In 
the  case  of  a  radome,  the  surfaces  have  finite  dimensions  creating  edges.  With  the  advent  of 
low  observable  large  RCS  contributors  have  been  reduced  to  such  an  extent  that  edge  effects 


1-3 


now  become  a  major  scattering  source.  Thus,  edge  effects  become  an  increasingly  important 
part  of  the  analysis  of  the  overall  scattering  of  objects. 

1.2  Problem  Statement 

In  his  dissertation,  Collins  [4]  proposed  an  extension  to  Munk’s  [11]  work  to  accommo¬ 
date  edges  by  using  the  concept  of  a  “physical  basis  function”  (PBF).  Collins  concluded  that 
at  some  distance  from  an  edge  of  a  sufficiently  long  finite  periodic  structure,  the  perturba¬ 
tions  in  the  amplitude  of  the  currents  induced  by  an  incident  field  would  approach  a  steady 
state.  At  this  “central”  location  where  a  steady  state  amplitude  is  attained,  Floquet  analysis 
can  be  used  to  determine  the  impedance  or  admittance  of  the  elements.  Then,  the  central 
portion  of  the  array  is  coupled  to  the  edge  element  basis  functions  in  a  moment  method 
solution. 

The  large  computational  expense  of  the  moment  method  in  calculating  the  scattered 
fields  of  complex  scatterers  and  inhomogeneous  materials  prompted  this  research.  FEM’s 
ability  to  compensate  for  the  MoM’s  shortfalls  in  the  areas  described  above  made  it  the  likely 
candidate  for  this  investigation.  Thus,  the  objective  of  this  research  is  to  devise  a  Finite 
Element  Method  formulation  for  modeling  a  semi-infinite  periodic  array  by  incorporating 
the  concept  of  a  physical  basis  function  as  a  model  for  array  elements  sufficiently  far  from 
edges. 

1.3  Scope 

The  focus  of  this  thesis  is  to  establish  the  validity  of  the  finite  element  method  to  solve 
truncated  periodic  structures.  Thus,  in  order  to  perform  a  proof-of- concept,  the  simplest 
cases  were  considered.  In  this  case  it  is  appropriate  to  confine  the  solution  domain  to  two- 
dimensional  space.  Confining  the  geometry  to  two  dimensions  makes  it  sufficient  to  use  the 
scalar  wave  or  Helmholtz  equations  without  worrying  about  spurious  modes^. 

^Spurious  modes  are  artifacts  of  numerical  calculations  that  often  appear  in  the  solution  when  using  the 
vector  Finite  Element  method.  They  have  no  physical  meaning  and  are  wrong  answers  [19]. 
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The  structure  used  to  develop  the  concepts  of  this  research  consists  of  perfect  electrical 
conducting  (PEC)  cylinders.  Dielectric  materials  axe  not  used  other  than  the  perfectly 
matched  anisotropic  (PMA)  absorbing  layers  [15]  on  the  boundary  of  the  solution  domain. 

PMA  layers  are  used  to  enclose  the  solution  domain  as  opposed  to  an  analytic  boundary 
condition  [2]  in  order  to  minimize  the  number  of  finite  elements  in  the  solution  domain. 
The  PMA  absorbing  boundary  turns  out  to  be  a  better  choice  for  implementing  a  periodic 
boundary  of  a  planar  semi-infinite  periodic  structure  since  it  allows  a  rectangular  boundary. 
Rectangular  boundaries  enclose  long  slender  objects  without  leaving  a  lot  of  free  space  in  the 
solution  domain.  In  contrast,  the  analytic  boundary  method  requires  a  circular  boundary. 

1.4  Thesis  Organization 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  2  provides  the  electromag¬ 
netics  background  for  the  Finite  Element  Method  applied  to  periodic  structures.  Chapter  3 
explains  the  methodology  used  in  completing  the  research.  The  results  of  the  simulations  and 
its  validation  are  shown  in  Chapter  4.  Conclusions  from  this  FEM  research,  including  pro¬ 
posed  areas  for  future  study  are  stated  in  Chapter  5.  Appendix  A  includes  all  the  field  plots 
from  the  geometries  of  interest.  They  include  contour  and  surface  plots.  In  Appendix  B,  the 
modifications  to  the  codes  formulated  by  Pelosi,  et  al.  [12]  are  included,  as  well  as  the  mesh 
construction  file  of  the  structure  used  to  test  the  hypothesis  in  this  thesis. 
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II.  Summary  of  Current  Knowledge 


2.1  Overview 

The  digital  computer  allows  the  computational  electromagnetics  (CEM)  practitioner 
to  add  the  Finite  Element  Method  to  his  portfolio  of  numerical  tools  because  of  its  greater 
potential  in  solving  more  complex,  real-world  engineering  problems. 

This  chapter  presents  the  summary  of  current  knowledge  leading  to  the  finite  element 
method  formulation  of  the  solution  to  the  scattering  of  a  semi-infinite  periodic  structure  in 
two  dimensions. 

2.2  Literature  Review 

2.2.1  Infinite  Periodic  Structures.  The  study  of  ESS  or  infinitely  periodic  struc¬ 
tures  begins  by  finding  a  method  of  reducing  the  entire  array  into  one  representative  element 
or  reference  cell.  Floquet’s  Theory  provides  the  methodology.  When  a  plane  wave  is  inci¬ 
dent  on  an  infinitely  periodic  array,  the  induced  fields  at  the  elements  in  the  array  are  also 
periodic.  In  other  words,  the  amplitude  of  the  fields  generated  by  every  cell  in  the  array  are 
identical  to  the  reference  cell  except  for  a  linear  phase  shift  factor  [11].  For  an  incident  plane 
wave  of  the  form 


E*(R)  = 

_  ^.^iQ-j0{sxX+Syy+Szz)  ^  (2.1) 

where:  ej  is  the  polarization  vector 
s  is  the  propagation  vector, 

the  total  electric  field  at  a  point  on  the  qrnl^  cell  (the  cell  in  the  row  and  column 
of  the  array)  can  be  expressed  as  a  function  of  the  field  at  the  equivalent  location  on  the 
reference  cell  at  the  origin  as  [11] 


E(r,„)  = 


(2.2) 


2-1 


or  the  magnetic  field  as 

M(r,^)  =  (2.3) 


where:  is  the  period  in  the  x  direction 

Dz  is  the  period  in  the  z  direction 
P  is  the  wave  number. 

For  an  infinite  array  of  Hertzian  dipoles  (of  differential  length,  dl)  as  shown  in  Fig¬ 
ure  2.1,  antenna  theory  gives  us  the  magnetic  vector  potential  at  an  arbitrary  observation 
point  a  distance  R  from  the  dipole  at  the  row  and  column  as  [11] 


dA 


qm 


=  P- 


Att  R 


qm 


(2.4) 


where  p  is  the  unit  vector  in  the  direction  of  the  current  elements.  Note  that  the  observation 
point  could  also  be  out  of  the  plane  where  the  dipoles  are  located. 
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Taking  all  the  contributions  from  every  element  in  the  array,  and  relating  the  Floquet 
currents  with  respect  to  the  current  in  the  reference  element  gives 
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The  electric  field  directly  obtained  from  Equation  (2.5)  is  a  very  slowly  converging 
series.  However,  using  the  Poisson  sum  formula 


Y,  (me,)  =  r  Y  fit  +  oT), 


(2.6) 


and  Fourier  transforming  the  electric  field  to  the  spectral  domain,  Munk  [11]  derived  the 
far-field  electric  field  of  an  array  of  Hertzian  dipoles  in  a  lossless  medium,  with  the  reference 
element  shifted  a  distance  R'  from  the  origin,  as 


nidi  °°  °° 

dE{x,y,z)  =  ^  E  - - [{r-p)f-p\. 


2DxDz  k=-oon=-oo 


(2.7) 


For  finite  length  dipole  elements,  the  electric  field  is 


E{x,y,z)  = 


00  00  g~i/?R*r± 


E  E 


2DxDz  ;t=_oo  n=-oo 


— [{r±  •  p)r±  -  p]  I  (2.8) 

jRef.Elem. 


where:  p  is  the  unit  vector  along  the  dipole’s  length 
[(r±  •  p)r±  —  p]  is  the  polarization  vector 
r±  =  ®(s^  +  1^-)  ±  Wy  +  z{s:,  +  I;) 

Vy  =  —  (5^  +  whose  root  is  either  positive  real 

or  negative  imaginary 

k,  n  are  the  new  row  and  column  indices  in  the  spectral  domain. 
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Having  found  the  electric  field  radiated  by  the  dipoles  in  the  array,  the  mutual  coupling 
(or  impedance  matrix)  between  them  can  be  calculated.  Since  the  excitation  vector  entries 
are  known  from  the  incident  field,  the  system’s  response  can  be  solved  [9,18]. 

2.S.2  Semi-Infinite  Periodic  Structures.  Semi-infinite  periodic  structures  present 
a  different  problem.  Although  the  available  literature  remains  relatively  sparse,  several 
approaches  have  been  published. 

Wasylkiwskyj  [22],  in  1973,  analyzed  the  mutual  coupling  effects  in  a  semi-infinite 
array  where  he  used  the  Weiner-Hoph  factorization  procedure  to  to  solve  an  infinite  order 
difference  equation  of  the  currents  formulated  at  the  antenna  ports. 

Cwik  and  Mittra  [6]  investigated  curved  and  truncated  strip  arrays.  Their  technique 
requires  two  piecewise  approximations.  The  first  approximation  assumes  the  induced  cur¬ 
rents  in  every  element  to  be  identical  as  if  belonging  to  an  infinite  array.  Then  they  replaced 
the  edge  element  currents  with  currents  found  from  the  solution  of  a  very  small  finite  array. 
Unfortunately,  this  method  does  not  perform  well  when  the  angle  of  incidence  is  near  grazing 
because  it  fails  to  take  into  account  the  increased  significance  of  the  coupling  between  the 
edge  elements  and  the  inner  elements. 

Collins’  approach  to  the  edge  problem  involves  a  modification  to  the  Poisson  Sum 
Formula  used  in  the  analysis  of  infinite  periodic  arrays.  The  one-sided  Poisson  Sum  Formula 
that  he  developed  specifically  for  semi-infinite  periodic  arrays  is  derived  by  defining  F{mu}o) 
in  Equation  (2.6) 


...  to  be  the  product  of  an  infinite  domain  continuous  function  G(lOo)  and  the 
Heaviside  unit  step  function  shifted  by  ^  to  capture  the  entire  m  =  0  term  [5] . 

Following  an  inverse  Fourier  transformation,  the  one-sided  Poisson  Sum  Formula  can  be 
expressed  as 


q=:Q  k=—oo 


'  g{t  ±  kT) 


+ 


27r  i-, 


9{'r)e 


{t±kT-T) 


t  ±  kT  —  T 


dr. 


(2.9) 
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Note  that  the  second  term  in  the  right  hand  side  is  a  principal  value  integral  and  can  be 
evaluated  using  residue  calculus. 

In  the  two  dimensional  case  where  the  test  geometry  was  an  array  of  infinitely  long 
thin  wires  which  are  periodically  spaced  extending  all  the  way  to  infinity  (see  Figure  2.2), 
Collins  showed  that  there  will  be  some  distance  amplitudes  will  be  approximately  Over  this 
section,  the  array  where  Munk’s  periodic  characteristics,  to  the  edge 

An  analog  in  the  Finite  Element  Methods  to  the  technique  just  described  is  what 
motivated  this  research. 

2.2.S  The  Finite  Element  Method  Model  of  a  Periodic  Structure.  Like  all  numerical 
computation  algorithms,  the  finite  element  method  requires  a  problem  domain  that  is  finite. 
For  the  FEM  however,  this  requirement  is  due  to  the  physical  discretization  of  the  domain. 
In  order  to  fully  enclose  a  structure  that  is  infinitely  long  but  periodic,  Floquet  analysis  is 
utilized.  Floquet  theory  serves  as  the  tool  for  collapsing  an  entire  array  into  a  representative 
cell,  similar  to  the  one  previously  described  in  the  method  of  moments. 

Implementation  of  the  FEM  on  periodic  structures  has  already  been  accomplished  by 
others.  Gedney,  et  al.  [8]  characterized  periodic  gratings,  while  McGrath  [10]  applied  the 
FEM  to  phased  array  antennas. 

Referring  to  Figure  2.3,  a  unit  cell  representing  a  periodic  element  of  an  infinite  array 
is  enclosed  by  an  imaginary  rectangular  boundary,  d£l.  The  top  and  bottom  boundaries, 
dQ,u  and  represent  radiation  boundaries.  The  two  side  walls,  30,1  and  OOr,  represent 
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periodic  boundaries.  In  order  to  account  for  the  interaction  between  cells,  McGrath  compares 
the  periodic  boundary  to  “open- circuit”  waveguide  side  walls  [10].  The  mesh  elements  in 
the  computation  region,  however,  must  be  symmetric.  That  is  the  left  and  right  halves  of 
the  cell  must  be  a  mirror  image  of  each  other.  This  further  implies  that  the  nodes  on  the 
left  boundary  of  the  cell  must  have  the  same  vertical  coordinates  as  the  nodes  on  the  right 
boundary.  In  effect,  the  procedure  amounts  to  a  “wrapping  around”  of  the  left  boundary 
into  the  right  boundary,  analogous  to  performing  a  circulaj  convolution  in  digital  signal 
processing. 

The  method  used  by  McGrath  combines  the  Finite  Element  Method  with  Floquet 
Modal  Expansion.  Hybrid  methods,  especially  the  ones  in  combination  with  integral  methods 
like  the  MoM,  essentially  model  the  exact  representation  of  the  electromagnetic  fields  outside 
the  domain  of  the  problem.  Like  in  the  “pure”  MoM  [9],  this  makes  a  point  on  the  FEM 
boundary,  dCl,  interact  with  all  other  points  on  the  boundary.  This  non-local  behavior 
essentially  creates  a  dense  coupling  matrix,  thereby  requiring  more  computational  resources. 
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The  other  type  of  FEM  boundary  problem  implements  an  absorbing  boundary  enclos¬ 
ing  the  domain.  With  this  method,  test  points  (or  nodes)  in  the  the  domain,  0,  interact  only 
with  local  or  neighboring  points.  This  implementation  produces  a  sparse  coupling  matrix,  a 
big  plus  when  computational  efficiency  is  a  big  factor  in  solving  a  problem.  This  is  one  area 
where  the  FEM  is  superior  to  the  MoM  or  other  integral  methods. 

2.2.4  Software  Implementation  for  a  Two-Dimensional,  Semi-Infinite  Periodic  Struc¬ 
ture.  McGrath  developed  a  three-dimensional  finite  element  code  called  Phased  Array 
Antenna  Analysis  (PARANA)  to  solve  infinite  periodic  structures  as  phased  array  antenna 
geometries.  This  code  is  very  versatile  and  can  handle  different  types  of  problems  such  as 
two-port  waveguide  devices,  array  radiation,  reflection  and  transmission  from  infinitely  pe¬ 
riodic  structures,  and  array  scattering.  However,  downgrading  PARANA  in  order  to  handle 
a  two-dimensional  case  would  be  counter-productive. 

Pelosi,  et  ah,  in  their  book  [12],  released  a  finite  element  code  specifically  for  two- 
dimensional  problems.  Therefore,  it  is  the  belief  of  this  author  that  it  would  be  simpler 
to  modify  this  code  to  present  a  proof-of- concept  problem  involving  two-dimensional,  semi¬ 
infinite  periodic  geometries. 

2.3  Development  Approach 

Figure  2.4  illustrates  the  geometry  of  interest  in  this  research.  It  involves  a  semi-infinite 
periodic  geometry  in  two-dimensional  space.  The  problem  domain  consists  of  cylinders 
enclosed  in  Perfectly  Matched  Anisotropic  (PMA)  absorbing  layers. 

The  analysis  of  a  semi-infinite  periodic  structure  begins  by  considering  Collins’  imple¬ 
mentation  of  a  physical  basis  function.  The  FEM  model  consists  of  an  edge  element  or  cell 
followed  by  a  finite  number  of  similar  cells  extending  to  the  right  of  the  edge  element.  The 
finite  number  of  cells,  together  with  enclosing  the  geometry  with  PMA  absorbing  layers, 
takes  care  of  the  problem  of  having  an  infinite  calculation  domain  or  an  infinite  number  of 
unknowns. 
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Figure  2.4  Five- Cylinder  Periodic  Array  Enclosed  in  PMA  and  a  Periodic  Boundary 

To  take  into  account  the  fact  that  the  structure  actually  extends  indefinitely  to  the 
right,  a  boundary  there  must  be  one  that  simulates  a  infinitely  extending  periodic  bound¬ 
ary.  McGrath  provides  a  hint  as  to  how  to  characterize  this  boundary.  The  “open-circuit” 
boundary  method  previously  described  provides  the  means  for  taking  into  account  the  other 
cylinders  to  the  right  of  the  boundary.  However,  instead  of  folding  it  over  onto  the  left-hand 
side  boundary  of  a  periodic  cell  within  the  geometry,  the  wrapping-over  occurs  with  a  cell  im¬ 
mediately  to  the  right  of  the  periodic  boundary,  a  periodic  cell  not  within  the  bounds  of  the 
problem  domain,  fi.  Periodicity  at  the  boundary  takes  into  consideration  the  contribution 
of  more  periodic  scatterers  outside  the  region  il. 

Using  an  absorbing  boundary  in  the  FEM  formulation  accounts  for  the  fact  that  the 
enclosed  region,  ft,  as  a  whole,  is  not  a  periodic  one  although  its  cells  are. 

2.4  Summary 

This  chapter  presents  the  chronological  development  leading  to  the  solution  of  the 
semi-infinite  array  problem  using  the  method  of  moments  and  some  hybrids  of  the  MoM. 
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This  chapter  also  presents  a  parallel  solution  to  the  infinite  array  problem  using  the  finite 
element  method.  The  next  chapter  details  the  FEM  approach  to  solving  the  semi-infinite 
array  problem. 
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III.  Methodology 


3.1  Introduction 

The  solution  of  a  boundary  value  problem  using  the  Finite  Element  Method  can  be  suc¬ 
cessfully  implemented  by  stepping  through  a  very  structured  procedure.  After  the  definition 
of  the  boundary  value  problem,  the  solution  can  be  solved  in  five  steps  [12,21],  namely: 

1.  Discretize  the  problem  domain  into  finite  elements. 

2.  Apply  the  functional  (Wave  or  Helmholtz  equation)  over  each  element  in  the  region. 

3.  Apply  the  boundary  conditions  and  collate  the  contribution  from  each  finite  element 
in  the  region  to  form  a  global  or  coupling  matrix. 

4.  Solve  the  system  of  equations. 

5.  Process  the  data  to  extract  the  pertinent  information  and  create  a  visual  representation 
of  the  quantities  of  interest. 

The  first  step,  preprocessing,  is  to  mesh  the  region  where  the  solution  is  desired.  This 
can  be  accomplished  by  using  graphic  design  or  Computer-Aided  Design  (CAD)  software. 
This  is  a  mandatory  step  for  all  numerical  solutions.  This  step  is  where  the  coordinates  of 
nodes  or  points  in  the  region  are  passed  on  to  the  numerical  algorithm. 

The  second,  third  and  fourth  steps  belong  to  the  finite  element  method  algorithm. 
This  is  where  the  differential  equation  is  mapped  into  a  discrete  system  so  that  it  can  be 
numerically  processed  in  a  digital  computer.  In  the  second  step,  the  integro-differential 
equation  is  discretized  into  a  set  of  linear  equations  so  that  it  can  be  converted  into  matrix 
format.  It  is  also  in  this  step  where  the  integro-differential  equation  is  enforced  on  every 
element  of  the  grid. 

The  third  step  is  required  if  a  boundary  condition  must  be  explicitly  enforced.  For 
example,  a  Dirichlet  boundary  condition  occurs  when  at  the  surface  of  a  Perfect  Electrical 
Conductor  (PEC),  the  tangential  electric  fields  are  zero.  This  type  of  condition  must  be 
explicitly  enforced  in  the  algorithm.  A  Neumann  boundary  condition,  on  the  other  hand. 
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specifies  the  normal  derivative  at  the  boundary.  For  an  incident  plane  wave  with  a  transverse 
magnetic  (TM)  polarization,  the  normal  derivative  of  the  magnetic  fields  on  the  surface  of  a 
PEC  are  zero. 

The  fourth  step  is  choosing  a  solving  mechanism  for  the  system.  Galerkin’s  method  is 
the  preferred  method  in  the  finite  element  method.  This  mechanism  converts  a  differential 
equation  into  a  set  of  linear  equations  for  numerical  computation.  It  essentially  comes  down 
to  solving  a  matrix  equation  of  the  form 

Mix]  =  [6],  (3.1) 

The  matrix  [A]  is  the  called  the  coupling  matrix  where  its  elements  depend  on  the  inter¬ 
relationships  of  the  nodes  in  the  grid.  The  vector  [a;]  is  the  set  of  unknown  parameters,  and 
the  vector  [h]  is  a  vector  of  forcing  or  excitation  functions.  Solvers  then  essentially  invert  the 
matrix  [A]  and  multiply  it  by  [6]  to  get  the  solution  vector  [a;]. 

The  final  step  is  the  post-processing  step  where  the  solution  vector  is  processed  and 
presented  in  a  visual  format  that  is  conducive  to  analysis. 

3.2  Theoretical  Background 

Throughout  this  document,  the  fields  indicated  herein  are  assumed  to  be  time-harmonic 
with  an  dependence  which  is  suppressed.  Also,  the  analysis  assumes  a  transverse  mag¬ 
netic  TMz  polarized  incident  plane  wave  illumination  so  that  the  electric  fields  are  the 
unknown  quantities.  If  the  magnetic  field  equations  are  desired,  the  duality  principle  can  be 
used  on  the  equations  presented  here. 

In  two-dimensional  space,  where  a  cylindrical  scatterer  is  illuminated  by  a  plane  wave  of 
transverse  magnetic  (TM)  polarization  as  illustrated  in  Figure  3.1,  the  fields  in  the  scatterer’s 
vicinity  can  be  obtained  by  the  two-dimensional  wave  equation 


Figure  3.1  Two-Dimensional  Scattering  Configuration  for  Cylindrical  Scatterers 

where:  V*  =  x|;  +  y^, 

ko  is  the  free  space  wave  number, 

tr  is  the  relative  permittivity  of  the  medium, 

pir  is  the  relative  permeability  of  the  medium,  and 

/(a;,  y)  is  the  forcing  function  or  external  excitation. 

Equation  (3.2)  is  considered  as  a  strong  form  of  the  wave  equation  because  the  differ¬ 
ential  operating  on  the  unknown  is  second  order.  In  order  to  apply  the  finite  element  method 
to  Equation  (3.2),  the  order  of  the  differential  operator  must  be  reduced  to  one.  This  results 
in  the  weak  form  of  the  wave  equation. 

If  we  let  a  function  approximate  the  unknown  function  the  equality  in  Equa¬ 
tion  (3.2)  no  longer  holds  but  instead  results  in  a  non-zero  residual, 

r{x,  y)  =  Vf  l—VtE,,{x,  y)\  +  k^CrEzix,  y)  -  f{x,  y).  (3.3) 

\f^r  ) 

The  ideal  value  for  E^,  is  one  that  reduces  the  residual  to  zero  over  the  entire  domain,  Q. 
However,  that  is  usually  not  numerically  achievable  all  over  the  points  in  the  domain.  Instead 
it  is  more  practical  (and  can  be  more  easily  implemented)  to  find  an  approximation  Ez  that 
would  reduce  the  residual  to  a  minimum  value  at  all  points  in  0.  The  minimization  process 
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can  be  done  by  weighting  the  residual  over  M  subdomians.  Geometrically,  the  residue  must 
have  no  projections  onto  each  of  the  M  weighting  functions.  Mathematically,  taking  the 
inner  product  of  the  weighting  functions  with  the  residue  gives  the  best  approximation  to 
the  solution.  Thus, 

j  j^W{x,y)r{x,y)d^l  =  0.  (3.4) 

Substituting  Equation  (3.3)  into  Equation  (3.4)  and  dropping  the  functional  notation  of 
dependency  in  x  and  y,  the  wave  equation  becomes 

/ X 

It  is  usually  the  case  in  Galerkin’s  method  that  the  weighting  functions  and  the  basis 
functions  are  the  same. 

Applying  the  vector  identity  (Green’s  theorem  or  integration  by  parts), 

V  •  (qA)  —  qV  ■  A  +  Vq  ■  A,  (3.6) 

on  the  first  term  on  the  left-hand  side  of  Equation  (3.5)  results  in 

VEVt  •  —  =  Vt  ■  f  —WVA]  -  —VtW  •  VtE,.  (3.7) 

fir  Xf^T  j  Mr 

Substituting  Equation  (3.7)  into  Equation  (3.5)  results  in 

/  /  Vt  •  l—WVtEA  -  —VtW  ■  VtE,  +  klcrWE,  -W  ds^O.  (3.8) 

J  JQ  [  \flr  J  Mr 

Application  of  the  Divergence  Theorem, 

/  /  Vi  •  l—WVtEA  ds=  i  —W{VtE,  •  h)dl,  (3.9) 

J  Ja  \Mr  /  Hr 
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on  the  first  term  of  the  left-hand  side  of  Equation  (3.8)  results  in  the  weak  form  the  the  wave 
equation. 


-—VtW  •  V<4  +  khrWE,  -  Wf 


ds+  i  —W(n  •  VtEMl  =  0, 

Jdi}  iXt 


(3.10) 


where  50  represents  a  closed  contour  enclosing  the  domain,  and  h  is  the  outward  pointing 
unit  normal  vector  on  the  contour. 

In  free  space  and  in  the  absence  of  excitation  sources,  the  weak  form  of  the  wave 
equation  reduces  to  the  weak  form  of  the  Helmholtz  equation 


-—VW  •  +  k'^trWE, 

flrp 


ds 


—W{h  ■  VE,)dl  =  0. 

917  jJ/f 


(3.11) 


3.2.1  Formulation  of  a  Circular  Absorbing  Boundary.  The  finite  element  method, 
as  its  namesake  implies,  requires  a  geometric  domain  that  is  finite.  However,  that  does  not 
mean  that  the  domain  of  the  problem  must  also  be  finite  [17].  This  important  distinction 
enables  the  FEM  to  solve  problems  with  unbounded  regions,  such  as  scattering  and  radiation 
problems.  Typically,  propagation  problems  are  solved  by  seeking  the  solution  to  the  wave 
equation  or  the  Helmholtz  equation. 

The  physical  interpretation  of  the  solution  to  the  wave  equation  or  the  Helmholtz  equa¬ 
tion  is  an  outward  traveling  wave,  where  the  fields  asymptotically  decay  to  zero  at  infinity. 
Applying  the  FEM  to  such  problems  already  imposes  difficulties  in  its  implementation.  The 
biggest  hurdle  is  that  the  domain  must  be  finite.  The  next  hurdle  is  how  can  an  infinite  do¬ 
main  problem  be  simulated  by  a  numerical  algorithm  that  requires  a  finite  solution  domain 
and  still  capture  the  characteristics  of  an  outwardly  expanding  solution. 

In  order  to  overcome  both  difficulties,  an  absorbing  boundary  condition  (ABC)  must  be 
specified.  Absorbing  boundaries  accomplish  the  task  of  absorbing  all  scattered  energy  from 
the  enclosed  scatterer  (similar  to  a  black  body)  and  effectively  preventing  reflected  energy 
back  into  the  interior  region  or  the  solution  domain.  It  also  acts  as  a  bounding  surface  for 
the  mesh  elements. 
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Figure  3.2  Cylindrical  Scatterer  Enclosed  in  a  Fictitious  Circular  Absorbing  Boundary 

Generally,  analytic  absorbing  boundaries  represent  the  Green’s  integral  at  the  bound¬ 
ary  with  an  a  priori  approximation.  For  instance,  Bayliss  and  Turkel  applied  the  a  pri¬ 
ori  knowledge  of  the  Hankel  function  approximation  for  cylindrical  wave  propagation.  In 
hybrid  methods,  however,  the  Green’s  integral  on  the  boundary  is  numerically  evaluated 
directly  [17].  This  distinction  explains  why  absorbing  boundary  conditions  preserve  the 
sparsity  of  the  coupling  matrix.  On  the  other  hand,  one  important  disadvantage  of  the  ab¬ 
sorbing  boundary  condition  is  that  its  accuracy  drops  off  when  the  a  priori  approximation 
is  no  longer  satisfied. 

The  development  of  absorbing  boundaries  in  FEM  is  generally  credited  to  Engquist 
and  Majda(see  [7])  who  introduced  the  basic  concepts  using  a  planar  boundary.  Bayliss  and 
Turkel  [2]  built  on  Engquist  and  Majda’s  ideas  and  proposed  a  circular  absorbing  bound¬ 
ary.  Bayliss  and  Turkel  developed  a  second  order  differential  operator  by  starting  from  the 
assymptotic  Hankel  expansion  approximation  of  an  outward  traveling  electromagnetic  field. 
Since  the  Hankel  functions  characterize  cylindrical  wave  propagation,  Bayliss  and  Turkel’s 
model  requires  a  circular  absorbing  boundary.  Figure  3.2  illustrates  such  implementation. 
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Incorporating  the  second-order  Bayliss  and  Turkel  ABC  model  into  the  weak  Helmholtz 
functional,  Equation  (3.11),  and  decomposing  the  total  field  into  its  scattered  and  incident 
field  components  results  in  (see  [12,13,21]  for  derivation) 


r  r  \  r  r  r  (iW  fiF^ 

IL  -  li 

where:  a{p)  =  ih  +  T,  -  ^ 

M  =  5(7%  • 

From  this  point  on,  the  V  operator  indicates  two-dimensional  differentiation  and  the  un¬ 
known  function  represents  the  approximating  function  introduced  in  the  previous  section 
without  the  tilde  notation. 

Discretizing  Equation  (3.12)  above  would  result  in  a  linear  equation  of  the  form  shown 
previously  in  Equation  (3.1).  Note  that  the  right  hand  side  of  Equation  (3.12)  is  associated 
with  the  incident  field,  which  from  a  scattering  problem  point  of  view  is  a  known  parameter. 

S.2.S  Formulation  of  a  Perfectly  Matched  Anisotropic  Absorbing  Boundary.  The 
circular  absorbing  boundary  is  a  major  development  in  the  implementation  of  the  FEM  to 
electromagnetic  scattering  analysis.  The  only  restriction  in  its  application  is  that  it  must  be 
strictly  circular  in  shape.  That  could  be  a  drawback  when  the  scatterer  being  enclosed  is 
elongated.  For  elongated  scatterers,  a  circular  boundary  requires  a  larger  radius,  creating  a 
larger  area  with  nothing  but  free  space.  Although  the  coupling  matrix  remains  sparse  due 
to  the  “nearest  neighbor”  or  local  nature  of  the  coupling  between  the  nodes  on  the  grid,  its 
size  would  require  computation  on  free  space  areas  which  are  not  of  interest. 

For  problems  regarding  slender  or  elongated  bodies,  such  as  aircraft,  a  rectangular 
boundary  is  a  logical  geometry.  A  computational  boundary  analogous  to  that  of  absorbers 
used  in  anechoic  chambers  was  formulated.  Berenger  [3]  first  proposed  the  idea  of  a  rect¬ 
angular  enclosure  for  an  absorbing  boundary  for  use  in  the  Finite  Difference,  Time-Domain 
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(FDTD)  numerical  method.  Berenger’s  technique,  however,  modifies  Maxwell’s  equations  to 
accommodate  FDTD’s  time  stepping  routines.  Sacks,  et  al.  [15]  then  modified  Berenger’s 
formulation  and  used  anisotropic  materials  in  order  to  implement  the  formulation  in  FEM 
without  the  modifications  to  Maxwell’s  equations. 

The  idea  is  to  enclose  the  scatterer  with  a  rectangular  boundary  located  in  the  near 
zone  similar  to  the  circular  boundary  (see  Figure  3.3).  The  artificial  absorbing  materials 
used  are  anisotropic  and  usually  implemented  with  layers  of  absorbers  above  a  conducting 
surface,  much  like  the  absorber  configurations  in  an  anechoic  chamber.  A  perfectly  matched 
anisotropic  (PMA)  absorber  features  perfect  absorption  of  any  arbitrarily  polarized  plane 
wave  of  infinite  frequency  bandwidth  and  any  angle  of  incidence.  Although  the  interface 
between  free  space  and  the  anisotropic  material  is  reflectionless,  the  conducting  surface  which 
backs  the  finite  set  of  PMA  layers  would  certainly  produce  a  reflection.  However,  careful 
selection  of  the  layer  thickness  and  e  and  ii  of  the  material  reduces  reflections  back  into  the 
computational  domain  to  the  point  where  they  don’t  perturb  the  solution  significantly. 
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Starting  from  Equation  (3.11),  the  fields  in  the  interior  region,  fi,  can  be  decomposed 
into  two  subregions,  the  free  space  region  and  the  PMA  region.  In  the  free  space  region  the 
weak  Helmholtz  equation  is 


IL 


— VVE  •  VEjn  -If  klerWE^jn  +  i  W{h  ■  VE,)dl  =  0, 
jj/fp  J 


(3.13) 


where:  U  —  VIa  denotes  the  free  space  computational  region 
CIa  is  the  computational  region  in  the  PMA  absorber. 

Similarly,  the  Helmholtz  equation  characterizing  the  fields  in  the  PMA  absorber  region  is 


J -^VW-VEM-  J  + W{h-VE,)dl  =  0,  (3.14) 


where  [jj]  and  [e]  are  tensors  characterizing  the  PMA  and  are  of  the  form  [15] 


[i^  =  f^o 


a  0  0 
0  &  0 
0  0c 


(3.15) 


and 


[e]  =  Co 


a  0  0 

0  b  0 

0  0c 


(3.16) 


The  diagonal  elements  a,  b,  and  c  in  the  permittivity  and  permeability  tensors  are  complex 
numbers.  For  perfect  impedance  matching  between  free  space  and  the  anisotropic  medium 


S  =  S  =  [A] 

fio 


a  0  0 
0  6  0 
0  0c 


(3.17) 


Referring  to  Figure  (3.3),  the  contour  integral  in  Equation  (3.14)  includes  two  bound¬ 
aries,  the  PEC  backing  and  the  interface  between  free  space  and  the  PMA.  The  contour 
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integral  at  the  PEC  vanishes  because  of  the  tangential  electric  fields  there  are  zero.  The 
other  part  of  the  contour  integral  in  Equation  (3.14)  is  the  same  as  the  contour  integral  of 
Equation  (3.13).  However,  since  their  surface  normals  are  opposite  in  direction,  they  cancel 
out.  Hence,  Equation  (3.13)  becomes 


—VW  •  - 

A  Mr 

—WtW  •  VtEidn  + 

fir 


1 1  k%wEida 

[  I  kltrWEidn, 
J  JQ-Qa 


(3.18) 


after  decomposing  the  total  field  into  its  scattered  and  incident  components.  Similarly,  the 
field  equation  in  the  PMA  region  is 


dWdEl 
^  dy  dy 


dQ,  + 


// 

J  JQa 


klcWEidO. 


b- — ^dn 

OX  ax 

V,W  ■  V,Eidn  +  J  klWEldn, 


(3.19) 


where  a,  b,  and  c  are  the  diagonal  elements  of  the  complex  permittivity  and  permeability 
tensors  characterizing  the  PMA  material.  Note  that  the  material  parameters  on  the  right 
hand  side  of  Equation  (3.19)  axe  set  to  unity  because  the  PMA  exists  only  for  the  scattered 
field  and  not  for  the  incident  field. 


3.3  Finite  Element  Grid  Generation  of  the  Problem 

Discretizing  or  meshing  the  problem  domain  is  usually  done  by  Computer-Aided  Draft¬ 
ing  (CAD)  packages.  The  software  package  that  was  utilized  [12]  was  not  sufficiently  au¬ 
tomated,  nor  was  it  enclosed  in  a  graphical  user  interface  (GUI).  The  problem  domain  is 
constructed  with  a  set  of  non-overlapping  quadrilateral  blocks.  A  block  is  defined  as  a  four¬ 
sided  object  in  two-dimensional  space  whose  sides  can  be  characterized  by  second-order 
curves  [12,16].  A  side  can  be  specified  by  a  set  of  two  end  points  and  an  intermediate  point. 
As  seen  in  Figure  3.4,  a  block  can  thus  be  defined  with  a  total  of  eight  points.  The  blocks 
are  further  subdivided  into  smaller  triangular  mesh  elements  in  the  mesh  generator. 
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Figure  3.4  Generic  Quadrilateral  Block  Description  [12] 


The  geometry  input  file  to  the  mesh  generator  requires  three  sections  to  fully  define  the 
domain  of  the  problem.  The  first  section  describes  the  set  of  non-overlapping  quadrilateral 
blocks,  their  node  numbering  sequence  based  on  a  global  node  numbering  system  for  the 
whole  geometry,  and  the  material  assignment  of  each  block.  Description  of  the  block  requires 
the  sequential  order  of  nodes  around  the  block  in  a  counterclockwise  direction.  The  second 
section  specifies  the  Cartesian  coordinates  of  each  node  in  the  region.  The  last  section 
describes  how  each  particular  block  is  to  be  subdivided  into  finer  sections.  The  geometry 
input  file  used  in  this  research  is  included  in  the  Appendix  B. 

The  output  of  the  mesh  generator  is  the  file  that  the  FEM  code  reads.  This  file  consists 
of  a  header  and  four  sections.  The  first  section  specifies  the  list  of  the  mesh  elements,  either 
triangular  or  quadrilateral,  the  nodes  that  make  up  the  element  (three  for  first-order  triangles 
and  four  for  first-order  quadrilaterals),  and  a  flag  indicating  which  material  the  mesh  elements 
belong  to.  The  second  section  specifies  the  nodes  and  their  Cartesian  coordinates  again,  as 
they  may  have  been  renumbered  by  the  optimization  routine  (i.e.  Delaunay  regularization) 
in  the  mesh  generator  [12].  The  third  and  fourth  sections  list  the  edges  and  edge  labels 
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Geometry 

Boundary 

Type 

Cylinder 

Diameter 

Incident  Angles 

Periodic 

lA 

0°  —  90°;  5°  increments 

PMA 

lA 

0°  —  90°;  5°  increments 

15-Cylinder 

PMA 

lA 

0°  —  90°;  5°  increments 

Table  3.1  Geometry  Descriptions  for  the  Analysis  of  a  Semi-Infinite  Periodic  Structure 

of  each  element  in  the  grid.  The  third  and  fourth  sections  are  particularly  useful  when 
using  vector  finite  elements  for  computing  three-dimensional  scattering.  However,  since  this 
research  is  confined  to  two-dimension  scattering,  edges  and  edge  labels  are  not  particularly 
used. 

Table  3.1  shows  the  geometries  of  interest  which  consist  of  three  sets,  each  made  up  of 
two-dimensional  cylindrical  scatterers.  The  first  set  is  a  semi-infinite  periodic  array  of  five 
cylinders  shown  in  Figure  2.4.  The  solution  domain  is  bounded  by  a  rectangular  enclosure 
consisting  of  Perfectly  Matched  Anisotropic  (PMA)  absorber  layers  on  three  sides  and  by  a 
periodic  boundary  on  one  side.  The  PMA  boundaries  enclose  the  domain  on  the  top  and 
bottom  sides  and  on  the  left  side  of  the  periodic  structure.  The  periodic  boundary  is  located 
on  the  right  side  of  the  structure.  The  cylinders  are  periodically  spaced  2 A  apart.  There  is 
also  a  five-cylinder  finite  array  that  is  fully  enclosed  by  a  rectangular  PMA  boundary.  Both 
geometries  have  cylinders  that  are  lA  in  diameter. 

The  other  geometry  set  is  a  finite  array  of  fifteen  cylinders,  shown  in  Figure  3.5.  The 
cylinders  are  also  periodically  spaced  2A  apart.  The  fifteen-cylinder  array  is  fully  enclosed 
by  a  rectangular  box  of  PMA  absorbing  layers. 

The  fields  that  are  scattered  by  the  5-cylinder  and  the  15-cylinder  finite  arrays  are 
used  for  comparison  with  the  fields  scattered  by  the  5-cylinder  periodic  array.  By  comparing 
the  fields  in  the  vicinity  of  the  edge  cylinders  with  the  fields  in  the  vicinity  of  the  cylinder 
next  to  the  periodic  boundary,  a  measure  of  the  effectiveness  of  the  implemented  periodic 
boundary  is  obtained. 

To  mesh  each  one  of  the  specified  geometries,  intermediate  grid  points  are  defined 
and  the  problem  domain  divided  into  blocks.  Intermediate  grid  points  are  usually  defined 
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Figure  3.5  Geometry  Description  for  15-Cylinder  Finite  Array  Fully  Enclosed  in  a  Rect¬ 
angular  PMA  Boundary 

on  interfaces  between  two  materials,  like  the  interface  between  the  PMA  and  the  free  space 
inside  the  problem  domain,  or  on  PEC  boundaries.  This  graphical  description  of  the  problem 
domain  is  the  required  input  to  the  mesh  generator. 

The  grid  or  mesh  generator  can  mesh  the  domain  with  elements  that  are  either  triangles 
or  quadrilaterals.  Triangles  are  usually  preferred  because  they  can  readily  approximate  any 
complex  shape.  The  degree  of  approximation  can  be  adjusted  by  increasing  or  decreasing  the 
number  of  triangles  in  the  mesh.  For  example,  to  accurately  render  an  arbitrarily  complex 
shaped  domain,  smaller  triangles,  and  therefore  more  of  them,  are  required.  More  triangles 
mean  more  nodes  (or  vertices  of  the  triangles)  which  translate  to  more  unknowns  required  to 
compute  the  solution  to  the  fields  in  the  domain.  Alternatively,  a  coarser  approximation  to 
the  same  domain  requires  fewer  triangles,  which  translates  to  a  smaller  number  of  unknowns. 
The  decision  to  increase  or  decrease  the  number  of  mesh  elements  depends  on  several  factors. 
These  are  the  computational  resources  available,  the  desired  accuracy  of  the  solution,  and 
the  time  it  takes  to  obtain  the  solution. 

In  the  finite  element  method  applied  to  electromagnetic  field  problems,  a  rule  of  thumb 
is  to  have  the  lengths  of  the  sides  of  triangles  about  ^  [16,21].  For  this  research,  the  domain 
is  meshed  with  triangle  side  lengths  of  |.  This  decision  is  based  on  consideration  of  the 
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Figure  3.6  Shape  Functions  and  Linear  Approximations  of  Fields  Inside  a  Triangular  Mesh 
Element  [12] 

time  it  takes  to  mesh  the  problem  domain,  especially  with  the  fifteen-cylinder  finite  periodic 
array.  The  other  factor  affecting  the  decision  to  use  |  is  to  minimize  the  number  of  nodes 
or  unknowns  in  the  problem  domain. 

S.4  Discretization  of  the  Weak  Form  of  the  Wave  Equation 

Discretizing  the  weak  Helmholtz  or  wave  equations  by  representing  them  as  a  linear 
system  of  equations  comes  after  meshing  the  problem  domain.  This  can  be  done  by  approxi¬ 
mating  the  fields  inside  each  mesh  element  with  a  set  of  linear  basis  functions,  also  commonly 
known  as  shape  functions.  For  first-order  tricingles,  the  approximating  function  is 

=  (3.20) 

i=l 

where  the  a^^^’s  are  the  basis  functions  or  shape  functions  for  each  of  the  nodes  of  the 
element  in  the  mesh  (see  Figure  3.6).  The  triangles  used  throughout  this  reaseach  are  first 
order,  meaning  they  have  a  total  of  three  nodes.  are  the  yet  unknown  coefficients  of 
the  basis  functions  and  they  represent  the  values  of  the  fields  at  each  node  of  the  mesh 
element. 
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Shape  functions  approximate  geometrical  patterns  with  polynomial  expressions.  For 
triangular  mesh  elements,  the  basis  or  shape  functions  are  usually  expressed  in  what  Silvester 
defines  as  simplex  coordinates.  Simplex  coordinates  transform  shape  functions  defined  in  the 
Cartesian  coordinate  system  to  a  coordinate  system  that  is  independent  of  position.  This 
is  important  since  the  triangular  elements  that  make  up  the  grid  do  not,  in  general,  have 
the  same  aspect  ratios  or  are  similarly  collocated  in  space.  For  example,  one  triangle  in  the 
mesh  may  have  a  different  orientation  compared  to  another  one  in  the  mesh.  The  simplex 
coordinate  transformation  equations  are  [16] 


at 
0^2 
03 

where  the  area  of  the  triangle  is: 

^  =  \[{^2  -  3;i)(j/3  -  y\)  -  (a;3  -  3:i)(i/2  -  yi)]- 
Further,  the  shape  functions  have  the  following  property 


(x2j/3  -  *32/2)  (2/2  -  2/3)  (*3  -  *2) 

1 

(xayi  -  Xiys)  (2/3  -  yi)  (*i  -  *3) 

X 

(*l2/2  -  *22/i)  (2/1  -  2/2)  (*2  -  *1) 

y 

I  1  iii  =  j 
I  0  Hi ^  j 


(3.21) 


(3.22) 


The  matrix  formulation  can  be  achieved  by  defining  a  weighting  or  testing  function  and 
talcing  its  inner  product  with  the  basis  functions.  The  most  preferred  method,  Galerkin’s 
method,  requires  the  weighting  function  to  be  equal  to  the  basis  functions. 

Note  that  the  only  reason  for  discretizing  an  integro-differential  equation  is  to  obtain 
N  linearly  independent  algebraic  equations  so  that  it  can  be  implemented  numerically  [14]. 
Galerkin’s  method  of  applying  weighting  functions  to  the  differential  equation,  and  making 
it  equal  to  the  basis  function  exactly  does  that  process. 
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Discretizing  the  Helmholtz  equations  in  the  the  PMA-enclosed  free  space  and  inside 
the  PMA  material  is  implemented  simply  by  substituting  the  approximating  functions  of 
Equation  (3.20)  into  Equations  (3.18)  and  (3.19).  Doing  so  results  in 

1(S(.)][£4.)]  +  1|SW]|£-C)]  -  (3.23) 

in  the  free  space  portion  of  the  computational  region  and 

inside  the  PMA  material.  The  right  hand  side  of  both  Equations  (3.23)  and  (3.24)  above  are 
associated  with  the  incident  fields  so  that  they  are  usually  known  quantities.  The  coefficient 
matrices  are 


[Si^)]=  I  f 

J  EaC®)  dx  dx 

(3.25) 

dy  a, 

(3.26) 

(3.27) 

The  fields  in  each  element  of  the  mesh  are  calculated  individually.  Generally,  inter¬ 
actions  or  coupling  between  nodes  occur  only  between  immediate  neighbors.  If  a  hybrid 
method  is  utilized,  however,  coupling  may  occur  between  nodes  other  than  its  immediate 
neighbor.  Thus,  when  assembling  the  coupling  or  global  matrix,  the  field  value  of  a  node 
from  a  triangle  is  summed  up  only  with  field  values  from  neighboring  nodes,  resulting  in  a 
sparse  matrix.  The  neighboring  nodes  may  either  belong  to  adjacent  triangles  or  from  the 
same  triangle. 

Part  of  the  assembly  process  for  the  coupling  matrix  is  the  incorporation  of  any  bound¬ 
ary  conditions.  A  Dirichlet  boundary  condition  exists  in  a  problem  domain  when  field  values 
are  prescribed  for  nodes  on  the  boundary.  For  instance  =  /(t),  where  f{t)  =  0  on  the 
surface  of  a  PEC  for  transverse  magnetic  (TM)  polarization.  Another  type  of  boundary 
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Figure  3.7  Two-Dimensional  Periodic  Structure  [12]. 

condition  is  the  Neumann  boundary  condition  where  =  g{t).  For  Transverse  Electric 
(TE)  polarization,  the  magnetic  field  on  the  surface  of  a  PEC  vanishes  so  that  =  0. 

3.5  Implementation  of  the  Periodic  Boundary 

The  main  objective  of  this  research  is  to  find  a  way  to  mimic  the  fields  generated  by  a 
semi-infinite  periodic  structure  illuminated  by  a  plane  wave  in  free  space  with  a  substitute 
geometry  that  is  finite.  The  technique  proposed  is  by  hybridizing  the  EEM  with  a  Floquet 
modal  expansion.  Briefly,  Floquet  analysis  takes  advantage  of  the  periodicity  in  an  array 
and  designates  a  periodic  element  or  unit  cell  in  an  array  to  represent  the  entire  array. 

For  a  two-dimensional,  infinitely  periodic  structure  (one  that  is  uniform  in  one  direction 
and  periodic  in  the  other)  as  shown  in  Figure  3.7,  the  total  fields  induced  by  an  incident 
T M2-polarized  plane  wave  on  the  array  can  be  expressed  as  [12] 

El{x  +  d,y)  =  El{x,y)e^^^^°^^’l‘^\  (3.28) 
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where  d  is  the  spatial  period,  k  is  the  propagation  constant  in  the  medium,  and  is  the 
incident  angle. 

Equation  (3.28)  says  that  the  amplitudes  of  the  fields  in  a  periodic  structure  are  es¬ 
sentially  the  same  as  that  of  some  reference  cell  located  at  an  arbitrary  reference  frame  (say 
the  origin  in  Cartesian  coordinates)  except  for  a  phase  factor. 

For  the  periodic  structure  in  Figure  3.7,  Equation  (3.11)  applies  to  the  region  inside  the 
fictitious  boundary.  The  contour  integral  for  a  rectangular  enclosure  comprises  four  segments, 
dilu,  dCli,  and  OCIl-  The  weak  Helmholtz  equation  inside  the  region  is  [10, 12] 


-—VW-VE,  +  khrWE, 

fir 

+  I  —W(h  •  VEMl  +  f  —W(h  ■  VEMl 

Jd^D  l^r  l^r 

+  /  —W(h  •  VEMl  =  0. 

JdQR  fir 


(3.29) 


The  top  boundary,  d^lu,  and  the  bottom  boundary,  dTo,  enforce  the  continuity  of  the 
tangential  fields  between  the  bounded  region  and  free  space  above  and  below  the  structure, 
respectively.  These  are  sometimes  referred  to  as  Floquet  Harmonics  [10].  The  left  boundary, 
dill,  and  right  boundary,  dCln,  enforce  the  periodic  boundary  and  are  related  to  each  other 
by  a  Floquet  phase  factor. 


/  —W(h-VEMy=  [ 

JdQji  fir  JdQi^ 


—W(h  ■ 

fir 


(3.30) 


In  order  to  easily  impose  the  periodic  boundary  conditions  as  described  previously, 
there  must  be  a  one-to-one  correspondence  of  nodes  on  each  side  of  the  periodic  boundaries. 
For  the  2-D  geometry  shown  in  Figure  3.7,  this  essentially  means  that  each  node  on  the  right 
boundary  must  have  the  same  y-coordinates  as  a  corresponding  node  on  the  left  boundary 
[10, 12].  This  procedure  is  viewed  as  what  McGrath  calls  “folding  over  the  boundaries  onto 
each  other  [10].” 

A  semi-infinite  periodic  structure  can  be  similarly  implemented.  The  top  and  bottom 
boundaries  are  terminated  by  PEC-backed  PMA  layers.  The  side  boundary  where  the  edge  is 
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located  is  also  terminated  with  PEC-backed  PMA  layers.  On  the  right  boundary,  however,  a 
periodic  boundary  is  implemented.  This  is  a  necessary  condition  to  account  for  the  cylinders 
extending  forever  to  the  right. 

As  can  be  seen  in  Figure  2.4,  the  folding  over  of  the  right  boundary  onto  the  left 
boundary  can  not  be  accomplished  on  this  particular  domain  because  the  array  elements 
are  not  representative  of  a  “unit  cell”  described  in  the  infinite  periodic  array  formulation. 
Periodicity  at  the  right  boundary  was  initially  thought  to  be  implemented  by  multiplying 
the  amplitudes  of  the  fields  there  by  the  Floquet  phase  factor,  ^ 
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Figure  3.8  Implementation  of  a  Single  Periodic  Boundary 

This  is  equivalent  to  attaching  a  fictitious  unit  cell  of  period  d,  similar  to  the  one 
previously  described  in  the  infinite  array  case,  adjacent  to  the  right  boundary  of  the  semi¬ 
infinite  structure.  Figure  3.8  illustrates  the  boundary  implementation.  The  left  boundary 
of  the  fictitious  unit  cell  is  related  to  its  right  boundary  by  Equation  (3.28).  To  enforce  the 
continuity  of  the  fields  across  this  boundary,  the  left  hand  boundary  of  the  fictitious  unit 
cell  is  made  equal  to  the  right  boundary  of  the  semi-infinite  geometry.  The  right  periodic 
boundary  is  then  implemented  in  Equations  (3.18)  and  (3.19)  where  they  become 
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//  —VtW-VtE^jn-ff  klerWEldn-^-  [  — mn  • 
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(3.32) 


It  is  important  to  note  that  the  periodic  boundary  conditions  in  Equations  (3.31) 
and  (3.32)  only  apply  to  the  boundary  and  not  to  the  entire  interior  region  in  the  domain, 
U.  The  discretization  of  Equation  (3.31)  and  Equation  (3.32)  can  be  converted  into  matrix 
format  by  using  Galerkin’s  method  previously  discussed. 

3.6  The  Finite  Element  Solution 

Using  Galerkin’s  method,  Equations  (3.31)  and  (3.32)  can  be  discretized.  The  resulting 
equations  are  similar  to  Equations  (3.23)  and  (3.24), 


Mr  Mr  ^ 


Jkdcos{<pi) 


[^(e) 


FreeSpace 


Jkdcos(<pi) 


(3.33) 


and 


Jkdcos{<l)i) 


(3.34) 
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Figure  3.9  Node  Interactions  in  an  Element  Matrix 


The  difference  is  when  the  node  being  considered  is  on  the  periodic  boundary.  The  elements 
of  matrices  [5^®)],  [5^®)],  and  that  interact  with  the  boundary  node  will  be  multiplied 

by  the  phase  factor  For  example,  in  Figure  3.9,  if  node  1  in  the  first  element 

belongs  to  the  periodic  boundary,  then  the  element  matrix  is 


[5(®=i)] 


^jkd  cos{4>i)g_^^  ^jkd  5l3 

522  ^23 

giMcos(^05'3,  532  533 


(3.35) 


The  matrices  [T^®^]  and  are  implemented  similarly  to  the  procedure  shown  for  [5^®)] 

above. 

All  the  processes  described  in  the  previous  section  are  incorporated  in  one  subroutine 
in  the  code  provided  by  Pelosi,  et  al.  [12].  The  rest  of  the  code  remained  unmodified.  Based 
on  the  flowchart  shown  in  Figure  3.10,  the  modified  subroutine  only  involves  the  branch 
where  the  absorbing  boundary  is  made  up  of  PMA  layers.  The  flowchart  for  the  modified 
subroutine  for  implementing  the  single  periodic  boundary  is  shown  in  Figure  3.11.  Code 
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modification  allowing  for  the  incorporation  of  a  semi-infinite  periodic  structure  using  the 
analytic  ABC  (circular  boundary)  is  not  implemented. 

The  FEM  software  utilizes  subroutines  from  the  freely  available  mathematical  libraries 
of  the  Linear  Algebra  Package  (LAPACK)  program  [1].  The  LAPACK  subroutines  are 
mainly  used  for  solving  the  linear  equation  shown  in  Equation  (3.1)  to  obtain  the  fields  in 
the  solution  domain. 

3.1  Post-processing  and  Data  Visualization 

The  output  of  the  FEM  consists  of  the  value  of  quantities  that  are  of  interest  in  each 
node  specified  on  the  grid.  In  this  research,  the  quantity  of  interest  is  the  electric  field  since 
we  were  interested  in  the  Transverse  Magnetic  (TM),  or  E-polarization  scattering  from  a 
two-dimensional  body.  However,  the  mesh  grid  is  usually  not  uniformly  spaced.  In  order  for 
graphical  interfaces  to  display  a  coherent  plot,  the  data  points  must  be  uniformly  spaced. 
Thus,  a  graphical  post-processor  performs  the  function  of  interpolating  between  the  data 
obtained  in  the  FEM  solution  and  extracting  uniformly  incremented  data  points.  The  post¬ 
processor  provided  in  the  software  package  is  utilized  without  any  modifications  to  the  code 
other  than  changing  the  parameters  to  accommodate  more  nodes  and  mesh  elements. 

A  Matlab  code  is  written  to  take  advantage  of  its  powerful  graphical  environment.  The 
output  of  the  post-processor  is  converted  to  Matlab  format  to  produce  the  contour  map  and 
the  3-dimensional  surface  plot  of  the  fields.  The  Matlab  code  is  shown  in  Appendix  B.3. 

3.8  Summary 

This  chapter  shows  the  formulation  used  to  modify  the  FEM  solution  for  regions  en¬ 
closed  by  PM  A  absorber  layers  to  include  semi-infinite  periodic  structures.  The  modifications 
involve  removal  of  a  PMA  boundary  on  the  opposite  side  of  the  edge  where  the  semi-infinite 
structure  is  “cut”  and  substituting  a  periodic  boundary  in  its  place.  The  periodic  bound¬ 
ary  is  formulated  from  the  hypothesis  that  the  periodic  fields  are  produced  from  a  spatially 
periodic  structure.  The  next  chapter  shows  the  results  of  the  formulations  described  in  this 
chapter. 
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Figure  3.10  Flowchart  for  the  FEM  Code. 
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IV.  Analysis  of  Results 

The  previous  chapter  outlined  the  theory  leading  to  the  calculation  of  the  fields  scattered  by 
a  semi-infinite  periodic  structure.  This  chapter  compares  the  FEM  results  from  a  periodic 
structure  with  two  finite  arrays  of  different  lengths.  The  objective  is  to  establish  the  validity 
of  substituting  a  periodic  boundary  in  place  of  the  truncated  part  of  a  semi-infinite  array. 
The  concept  of  the  physical  basis  function  is  used  to  verify  this. 


4.1  Comparison  of  a  Fifteen- Cylinder  Finite  Periodic  Array  With  A  Five-Cylinder  Finite 

Periodic  Array  at  90°  Incidence. 

The  total  field  magnitude  map  generated  by  a  plane  wave  illuminating  the  15-cylinder 
array  at  90°  incidence  is  shown  in  Figure  4.1.  Looking  at  the  contour  plot,  90°  incidence 
is  coming  from  the  top,  in  the  plane  of  the  paper.  Visually  inspecting  the  contour  map, 
the  uniformity  of  the  fields  at  the  central  portion  of  the  array  is  evident.  It  also  shows  the 
perturbations  at  the  edges.  Further  inward  from  the  edges,  the  fields  gradually  approach  a 
uniform  pattern  in  the  central  portion  of  the  array.  The  bottom  plot  of  Figure  4.1  shows  a 
surface  plot  of  the  fields. 

To  show  the  PBF  concept,  we  start  by  taking  the  fields  in  the  vicinity  of  the  central 
cylinder  (cylinder  8)  and  subtracting  the  fields  from  the  adjacent  cylinders  (7  or  9,  cylinder 
1  is  the  left  edge).  The  difference  in  field  magnitudes  between  cylinders  8  and  7  is  shown 
in  Figure  4.2,  and  the  difference  in  field  magnitudes  between  cylinders  8  and  9  is  shown  in 
Figure  4.3.  The  difference  plots  show  that  the  field  values  vary  by  no  more  than  1%  of  the 
field  magnitudes  in  the  vicinity  of  cylinder  This  comparison  is  done  by  taking  the  highest 
magnitude  of  the  field  in  the  vicinity  of  cylinder  8  and  the  difference  in  field  magnitudes 
between  cylinder  8  and  the  the  other  cylinders, 


max[Cyls  -  CyU 
max[Cyl^ 


X  100%. 


(4.1) 
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Outlying 

Periods 

Max.  Field 
Difference 

Percent 

Difference 

Cylinder  7 

0.021 

1.0 

Cylinder  6 

0.043 

2.0 

Cylinder  5 

0.070 

3.3 

Cylinder  4 

0.093 

4.4 

Cylinder  3 

0.130 

6.2 

Table  4.1  Percent  Difference  in  Field  Magnitudes  Between  the  Central  Period  of  the  15- 
Cylinder  Finite  Array  and  its  Outlying  Periods 

Table  4.1  shows  the  percentage  differences  between  the  central  cylinder  (cylinder  8) 
and  the  outlying  cylinders.  For  an  arbitrary  cutoff  of  5%  difference,  it  is  safe  to  say  that 
the  fields  after  cylinder  3  can  be  defined  as  the  “steady  state”  portion  of  the  fields  in  the 
array.  Therefore,  in  order  to  implement  a  semi-  infinite  cylinder,  the  truncation  point  or 
the  placement  of  the  periodic  boundary  must  occur  after  4  periods.  To  further  provide  a 
margin  for  numerical  error  as  well  as  a  “convenient”  round  number  ratio  between  the  number 
of  the  semi-infinite  array  is  modeled  with  an  additional  of  5  periods  before  the  truncation 
point.  Figures  4.4  and  4.5  show  the  differences  in  field  values  between  cylinders  8  and  6,  and 
cylinders  8  and  5  respectively. 
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Figure  4.2  Difference  in  Field  Values  Between  Cylinders  8  and  7  of  the  Fifteen-Cylinder 
Finite  Array  at  90°  Incidence 

Figure  4.6  shows  the  fields  from  a  five-cylinder  finite  array.  It  is  obvious  that  the  fields 
are  also  nearly  symmetric  with  respect  to  the  center  cylinder.  It  should  be  obvious  that 
subtracting  the  fields  of  a  five-cylinder  array  from  the  fields  of  a  fifteen- cylinder  array,  with 
the  left  edges  aligned,  would  result  in  a  large  variation  in  the  difference  plots  at  the  fifth 
cylinder.  Figure  4.7  shows  the  difference  in  field  values  for  the  first  five  cylinders  of  the 
15-cylinder  finite  array  and  the  5-cylinder  finite  array.  Figure  4.8  shows  the  difference  in 
field  values  at  the  fifth  cylinder.  Note  the  scale  difference  from  the  previous  plots. 

4.2  Comparison  of  a  Fifteen-Cylinder  Finite  Periodic  Array  With  A  Five-Cylinder  Semi- 

Infinite  Periodic  Array  at  90°  Incidence. 

The  purpose  of  implementating  a  periodic  boundary  on  the  truncated  part  of  a  periodic 
array  is  to  simulate  a  physical  continuity  of  periodic  elements.  It  should  also  eliminate  the 
edge  diffraction  that  would  otherwise  occur  at  the  discontinuity.  As  seen  from  the  previous 
section,  the  right  boundary  of  the  5-cylinder  array  showed  huge  differences  in  the  fields. 
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Figure  4.3  Diiference  in  Field  Values  Between  Cylinders  8  and  9  of  the  Fifteen- Cylinder 
Finite  Array  at  90°  Incidence 


Figure  4.4  Difference  in  Field  Values  Between  Cylinders  8  and  6  of  the  Fifteen-Cylinder 
Finite  Array  at  90°  Incidence 
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Figure  4.7 


Figure  4.8 
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Contour  Plot  and  Difference  in  Field  Values  Between  the  First  Five  Cylinders 
of  the  15- Cylinder  and  5- Cylinder  Finite  Arrays  at  90°  Incidence 


Cylinder  5  Field  Differences  Between  the  5-Cylinder  Finite  Array  and  the  15- 
Cylinder  at  90°  Incidence 
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Wavelengths 

Figure  4.9  DiflFerence  in  Cylinder  5  Field  Values  for  Both  15-Cylinder  Finite  and  5-Cylinder 
Semi-Infinite  Arrays  at  90°  Incidence. 

Taking  the  5-cylinder  finite  array  and  replacing  the  right  edge  with  a  periodic  boundary 
converts  the  geometry  into  a  semi-infinite  geometry  as  far  as  FEM  is  concerned. 

Figure  4.9  shows  the  total  fields  of  the  semi-infinite  array.  Visually  inspecting  the 
fields,  one  can  see  the  significant  change  in  the  field  behavior  near  the  periodic  boundary. 
Figure  4.11  shows  the  difference  in  field  values  between  cylinder  5  of  the  15-cylinder  fi¬ 
nite  array  and  cylinder  5  of  the  simulated  semi-infinite  array.  Comparing  Figure  4.11  with 
Figure  4.8,  the  periodic  boundary  performs  effectively. 

4-3  Other  Incidence  Angles 

The  results  presented  thus  far  involved  90°  incidence.  For  off-normal  incidence,  the 
periodic  boundary  does  not  perform  cis  well  as  at  normal  incidence.  This  is  evident  in  the 
field  plots  of  the  semi-infinite  array.  Taking  the  75°  incidence  as  an  example,  the  contour 
and  field  map  from  the  semi-infinite  array  is  shown  in  Figure  4.12.  The  field  and  contour 
map  from  the  finite  array  is  shown  in  Figure  4.13. 

Looking  at  the  field  patterns  in  the  vicinity  of  the  cylinder  next  to  the  periodic  bound¬ 
ary  (right  edge)  in  Figure  4.12,  and  comparing  that  with  the  field  pattern  one  can  see  the 
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Figure  4.10  Difference  in  Field  Values  for  Both  15-Cylinder  Finite  and  5-Cylinder  Semi- 
Infinite  Arrays  at  90°  Incidence. 


Figure  4.11  Difference  in  Cylinder  5  Field  Values  for  Both  15-Cylinder  Finite  and  5- 
Cylinder  Semi-Infinite  Arrays  at  90°  Incidence. 
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Wavelengths 


Figure  4.12  Contour  Plot  and  Total  Field  Map  of  a  5-Cylinder  Semi-Finite  Array  at  75° 
Incidence 

discrepancies  in  the  implemented  periodic  boundary.  Comparing  fifteen- cylinder  array  shown 
in  Figure  4.1  with  the  five-cylinder  can  immediately  recognize  the  disturbance  in  periodic 
boundary.  The  fields  do  not  with  a  periodic  array. 

The  difference  in  the  field  values  between  the  first  five  periods  of  the  fifteen-cylinder 
finite  array  and  semi-infinite  array  is  shown  in  Figure  4.14.  The  difference  in  the  field 
values  between  the  fifteen- cylinder  finite  array  and  the  five-cylinder  finite  array  is  shown  in 
Figure  4.15.  is  apparent  that  the  periodic  boundary  does  not  perform  as  expected. 

Several  factors  may  contribute  to  the  failure  of  the  periodic  most  likely  explanation  is 
that  the  boundary  is  the  solution  domain.  The  phase  factor,  at  the  periodic  boundary  would 
At  normal  incidence  and  therefore  While  at  other  incident  angles,  reflection  at  the  boundary 
occurs  and  energy  would  be  added  into  the  solution. 

4.4  Alternative  Implementation  of  the  Periodic  Boundary 

The  periodic  boundary  implementation  presented  here  suggests  that  the  performance 
of  an  “open-circuit”  boundary  is  not  sufficient  in  fully  characterizing  the  fields  in  a  semi¬ 
infinite  periodic  array.  The  next  logical  step,  given  enough  time,  would  be  to  implement  a 
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Figure  4.14 


Contour  Plot  and  Field  Difference  Plot  for  the  First  Five  Periods  of  the  15- 
Cylinder  and  5-Cylinder  Finite  Arrays  at  75°  Incidence 
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Figure  4.15  Contour  Plot  and  Field  Difference  Plot  for  the  First  Five  Periods  of  the  15- 
Cylinder  Finite  Array  and  the  Semi-Finite  Array  at  75°  Incidence 

“folding  over”  of  the  periodic  boundary  with  the  opposite  side  of  the  unit  cell  containing 
cylinder  5  (the  cylinder  next  to  the  periodic  boundary).  Referring  to  Figure  4.16,  the  field 
values  at  the  nodes  in  the  periodic  boundary  are  related  to  the  field  values  a  distance  of  one 
period  from  the  boundary. 
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Figure  4.16 


Implementation  of  a  Single  Periodic  Boundary  by  Relating  the  Fields  in  the 
Interior  One  Periodic  From  the  Periodic  Boundary 
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V.  Conclusion  and  Recommendations 

5.1  Conclusion 

The  objective  of  this  thesis  is  to  provide  a  framework  for  developing  a  Finite  Element 
Method  solution  to  semi-infinite  periodic  arrays  by  utilizing  the  concept  of  a  physical  basis 
function.  It  is  known  beforehand  that  the  amplitude  of  the  currents  induced  on  the  elements 
far  enough  from  the  edge  of  a  semi-periodic  array  are  constant.  Based  on  that  pretext,  a 
periodic  boundary  using  a  Floquet  expansion  is  implemented  at  a  sufficient  distance  from 
the  edge  of  a  semi-infinite  periodic  array.  However,  as  seen  in  the  previous  chapter,  the 
implementation  was  only  effective  at  normal  incidence. 

5.2  Recommendations  for  Further  Study 

A  next  logical  step  is  to  implement  the  periodic  boundary  by  relating  the  fields  there 
with  the  fields  one  interior  period  from  the  boundary  (the  truncation  point).  This  amounts 
to  performing  the  same  procedures  as  in  the  implementation  of  an  infinitely  periodic  array. 
However,  the  interior  “quasi-boundary”  exhibits  coupling  with  other  interior  nodes  outside 
the  unit  cell  (see  Figure  4.16). 

Another  approach  might  be  implemented  by  using  the  vector  finite  element  method 
(VFEM)  on  the  “open-circuit”  boundary  implementation.  VFEM  is  primarily  implemented 
to  fix  the  problem  of  spurious  modes,  but  could  be  used  here  since  the  periodic  boundary 
may  be  contributing  anomalous  solutions  at  off-normal  incidence. 

Although  the  length  of  the  periodic  array  used  is  sufficiently  long,  the  central  portion  of 
the  array  where  the  field  amplitudes  are  “constant”  is  narrow.  A  longer  steady  state  region 
would  ensure  that  the  edge  effects  are  fully  “damped”  by  the  time  it  reaches  the  central 
portion  of  the  array.  A  longer  array  would  mean  more  nodes  and  elements  in  the  finite 
element  solution  domain.  Meshing  usually  takes  a  long  time,  and  larger  problems  would 
mean  longer  mesh  times.  Thus,  a  more  efficient  or  optimized  mesh  generator  needs  to  be 
developed. 
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Usually  in  scattering  and  radiation  problems,  the  RCS  is  one  quantity  of  interest. 
RCS  calculations  can  be  done  by  taking  a  contour  integral  around  the  solution  region  and 
projected  into  the  far  field. 

Finally,  once  the  problem  of  calculating  off-normal  incidence  is  fixed,  several  applica¬ 
tions  can  be  implemented.  Since  the  FEM  accommodates  dielectrics  handily,  application  of 
material  treatments  for  metallic  edges  can  be  studied. 
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Appendix  A. 


Field  Plots  for  the  Fifteen- Cylinder  Finite  Periodic  Structure 
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Figure  A. 8  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1  —  A  Diameter,  15-Cylinder  Finite  Array  at 

35°  Incidence. 
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A. 2  Field  Plots  for  the  Five-Cylinder  Semi-Infinite  Periodic  Structure 
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Figure  A. 20  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  0°  Incidence. 
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Figure  A. 21  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  5°  Incidence. 
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Figure  A. 24  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of 
Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  20°  Incidence. 
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Figure  A.25  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of 
Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  25°  Incidence. 
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Figure  A. 28  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  40°  Incidence. 
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Figure  A. 29  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  45°  Incidence. 
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Figure  A. 30  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  50°  Incidence. 
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Figure  A. 31  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  55°  Incidence. 
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Figure  A. 32  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  60°  Incidence. 
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Figure  A. 33  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  65°  Incidence. 
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Figure  A. 34  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  70°  Incidence. 
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Figure  A. 35  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  75°  Incidence. 
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Figure  A. 38  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1  —  A 
Diameter,  5-Cylinder  Semi-Infinite  Periodic  Array  at  90°  Incidence. 
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A. 3  Field  Plots  for  the  Five- Cylinder  Finite  Periodic  Structure 
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Figure  A. 39  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5- Cylinder  Finite  Periodic  Array  at  0°  Incidence. 
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Figure  A.40  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5- Cylinder  Finite  Periodic  Array  at  5°  Incidence. 
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Figure  A. 41  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Finite  Periodic  Array  at  10°  Incidence. 
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Figure  A. 42  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5-Cylinder  Finite  Periodic  Array  at  15°  Incidence. 
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Figure  A.43  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Finite  Periodic  Array  at  20°  Incidence. 
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Figure  A.44  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Finite  Periodic  Array  at  25°  Incidence. 


Wavelengths 


klap  and  Three-Dimensional  Map  of  t 
5-Cylinder  Finite  Periodic  Array  at  3C 


Wavelengths 


Figure  A.47  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of 
Diameter,  5-Cylinder  Finite  Periodic  Array  at  40°  Incidence. 
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Figure  A. 48  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Finite  Periodic  Array  at  45°  Incidence. 
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Figure  A. 53  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of 

Diameter,  5-Cylinder  Finite  Periodic  Array  at  70°  Incidence. 
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Figure  A. 54  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Finite  Periodic  Array  at  75°  Incidence. 
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Figure  A. 55  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 
Diameter,  5-Cylinder  Finite  Periodic  Array  at  80°  Incidence. 
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Figure  A. 56  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1 

Diameter,  5- Cylinder  Finite  Periodic  Array  at  85°  Incidence. 
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Figure  A. 57  Contour  Map  and  Three-Dimensional  Map  of  the  Total  Fields  of  a  1  —  A 

Diameter,  5-Cylinder  Finite  Periodic  Array  at  90°  Incidence. 
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Appendix  B. 

B.l  Finite  Element  Method  Software  Modification 


subroutine  EMBEDB  (UPLO,MAXNODE,MAXNXE,MAXELE,MAXDIEL,MAXBAND, 

&  ELE , xy , NLAB , ELAB , LPMAX ,LPMAY , LPMAXY , APMA , LDIR , 

&  DIRC , COEFFP , COEFFQ , nele ,nnode , A , B , ku , 

&  kl,kd,phi,Ndie) 

C 

C  EMBEDS  LOCAL  MATRIX  SL  OF  THE  ELEMENT  IE  INTO  THE  GLOBAL  MATRIX 
C  A  OR  INTO  THE  RIGHT  HAND  SIDE,  AS  APPROPRIATE.  THIS  EXPLOITS 
C  PMA  TYPE  ABC. 

C  DUMMY  ARGUMENTS  ARE  COMMENTED  IN  THE  CALLING  MODULE 
C  QUICK.FEM  (C)  1997  PELOSI  -  COCCIOLI  -  SELLERI 

C  This  code  allows  the  inclusion  of  a  periodic  boundary  on  one  side 
C  and  PMA  layers  on  the  other  three  sides.  This  enables  computation 
C  of  a  HALF-INFINITE  periodic  cylinders. 

C 

C  Modified  to  accommodate  the  Right  Hand  Side  Boundary  not  enclosed 

C  by  PMA.  Also  hard  coded  to  run  on  a  specific  geometry: 

C  — >  f ivepecyl.ren.f em 

C 

C  Modification  only  occurs  in  the  subroutine:  EMBEDB 

C 

C 

C  PERRY  N.  VILLANUEVA 

C  Air  Force  Institute  of  Technology 

C  6  January  1999 

C 

C  [IN] 

IMPLICIT  NONE 

CHARACTER*!  UPLO  ! ’U’  Stores  only  Upper  Triangle 

!’L’  Stores  only  Lower  Triangle 
!’T’  Stores  all  matrix 

integer  MAXNODE , HAXNXE , MAXELE , MAXDIEL , MAXBAND , ELE (0 : MAXNXE ,MAXELE) 
integer  NLAB (MAXNODE) , ELAB (MAXELE) 

integer  LPMAX, LPMAY, LPMAXY  !  Label  of  elements  on  the  PMA 
integer  LDIR  !  Label  of  nodes  on  Dirichlet  boundary 

integer  DIRC  !  Coefficient  of  scattered  field 

!  at  Dirichlet  boundeiry  (Uscat=DIRC*Uinc) 
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complex  CoeffP(MAXDIEL) ,CoeffQ(MAXDIEL) ,APMA,AA,BB,CC 

integer  Nele , NNode , Ndie , ku , kl ,kd , ie 

real  xy(2,MAXN0DE) ,SXE(8,8) ,SYE(8,8) ,TE(8,8) ,phi 

C  CIN/OUT] 

complex  A(3*MAXBAND+1,MAXN0DE)  !  Left  hand  side  matrix 

complex  B(HAXNODE)  !  Right  Hand  side  vector 

C  [LOCAL] 

integer  I,J  !  indices 

integer  IR0W,IC0L  !  indexes  to  entries  global  FEM  Matrix 

complex  caux 
real  PI,K0,K02 

C  [EXTERNAL  FUNCTIONS] 

complex  CBGET 
complex  EINC 

parameter  (PI  =  3.141592653589793238) 

KO  =  2.  *  PI 

K02=  KO  *  KO 

do  ie=l,nele 

call  ELEAMAT (MAXNODE , MAXNXE , MAXELE , XY , ELE , IE , SXE , SYE , 8 , TE , 8) 

DO  I  =  1,3 

IRON  =  ELEd.IE) 

IF  (NLAB(IROW)  .EQ.  LDIR)  THEN 


C  -  ENFORCE  DIRICHLET  BOUNDARY  CONDITIONS 

call  CBPUT((1.,0.), 

&  A,IR0W,IR0W,ku,kl,kd,3*MAXBAND+l,NN0DE,’U') 

if  (ELAB(ie) .ne.LPMAX.and.ELAB(ie) .ne .LPMAY. and. 

&  ELAB(ie) .ne.LPMAXY)  then 

B(IROW)  =  DIRC  *  EINC(Xy(l, IRON) ,xY(2,IR0W), PHI) 
endif 


ELSE 

DO  J  =  1,3 

ICOL  =  ELE (J, IE) 

IF  (ELAB(IE).le.NDIE)  THEN 


C  - AUGMENT  THE  GLOBAL  MATRIX  A 

C 

C  - >N0DES  ARE  IN  AIR  (INTERIOR) 

caux  =  CBGET(A,irow,icol,ku,kl,kd,3*MAXBAND+l, 
&  NNODE, 'U’) 
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if  (  (icol.eq.l895) 
(icol.eq. 1911) 
(icol.eq. 1926) 
(icol.eq. 1952) 
(icol.eq. 1981) 
(icol.eq. 2010) 
(icol.eq. 2042) 
(icol.eq. 2073) 
(icol.eq. 2078) 


.or.  (icol.eq. 1897) 
.or.  (icol.eq. 1913) 
.or.  (icol.eq. 1939) 
.or.  (icol.eq. 1966) 
.or.  (icol.eq. 1995) 
.or.  (icol.eq. 2026) 
.or.  (icol.eq. 2057) 
.or.  (icol.eq. 2077) 
)  then 


.or. 

.or. 

.or. 

.or. 

.or. 

.or. 

.or. 

.or. 


caux  =  caux  +  Coeffp(ELAB(IE)) 

*  (SXE(i,j)  +  SYE(i,j)) 

*  exp(-(0.,  I.)*k0*cos(phi)) 

-  k02  *  C0EFFq(ELAB(IE))  *  TE(i,j) 

else 

caux  =  caux  +  Coeffp(ELAB(IE)) 

*  (SXE(i,j)  +  SYE(i,j)) 

-  k02  *  C0EFFq(ELAB(IE))  *  TE(i,j) 

endif 

call  CBPUT(caux, 

A , IRO W , ICOL , ku , kl , kd , 3*MAXBAND+ 1 , NNODE , ’ U ’ ) 


AUGMENT  THE  RIGHT  HAND  SIDE 


if  (  (icol.eq. 1895) 
(icol.eq. 1911) 
(icol.eq. 1926) 
(icol.eq. 1952) 
(icol.eq. 1981) 
(icol.eq. 2010) 
(icol.eq. 2042) 
(icol.eq. 2073) 
(icol.eq. 2078) 


.or.  (icol.eq. 1897) 
.or.  (icol.eq. 1913) 
.or.  (icol.eq. 1939) 
.or.  (icol.eq. 1966) 
.or.  (icol.eq. 1995) 
.or.  (icol.eq. 2026) 
.or.  (icol.eq. 2057) 
.or.  (icol.eq. 2077) 
)  then 


.or. 

.or. 

.or. 

.or. 

.or. 

.or. 

.or. 

.or. 


B(IROW)  =  B(IROW)  -  (coeffp(ELAB(IE))* 

(SXE(i,j)  +  SYE(i,j)) 

*exp(-(0.,  I.)*k0*cos(phi)) 

-  K02*COEFFq(ELAB(IE))  * 

TE(I,J))  *  EINC(Xy(l, ICOL), xY(2, ICOL), PHI) 

else 

B(IROW)  =  B(IROW)  -  (coeffp(ELAB(IE))* 

(SXE(i,j)  +  SYE(i,j))  -  K02* 
C0EFFq(ELAB(IE))  * 

TE(I,J))  *  EINC(Xy(l, ICOL), xY(2, ICOL), PHI) 


endif 

ELSE 

—AUGMENT  THE  GLOBAL  MATRIX  A 

- >  NODES  ARE  IN  PMA 

if(ELAB(IE).EQ.LPMAX)  then 
AA=APMA 
BB=1./APMA 
CC=APMA 

else  if(ELAB(IE) .EQ.LPMAY)  then 
AA=1./APMA 
BB=APMA 
CC=APMA 
else 
AA=1. 

BB=1. 

CC=APMA*APMA 

endif 

caux  =  CBGET(A,irow,icol,ku,kl,kd,3*MAXBAND+l, 
NNODE.’U’) 

if  (  (icol.eq. 1880)  .or.  (icol.eq. 1894)  .or. 

(icol,eq.2082)  .or.  (icol.eq. 2085)  )  then 
caux  =  caux  + 

(SXE(i,j)/BB  +  SYE(i,j)/AA) 

*  exp(-(0.,  I.)*k0*cos(phi)) 

-  k02  *  CC  *  TE(i,j) 

else 

caux  =  caux  + 

SXE(i,j)/BB  + 

SYE(i,j)/AA  + 

-  k02  *  CC  *  TE(i,j) 

endif 

call  CBPUT(caux, 

A , IRO W , ICOL , ku , kl , kd , 3*MAXBAND+ 1 , NNODE , ’ U ' ) 


if  (nlab(ele(i,ie)) .eq.O)  then 

if  (  (icol.eq. 1880)  .or.  (icol.eq. 1894)  .or. 
(icol.eq. 2082)  .or.  (icol.eq. 2085)  )  then 

B(IROW)  =  B(IROW)  - 

((SXE(i,j)  +  SYE(i.j)) 


RETURN 

END 
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B.2  Five-Cylinder  Periodic  Array  Mesh  Description  File 


f ivepecyl_open . f em 
133  33 


1 

3 

1 

2 

3 

15 

23 

22 

21 

14 

55 

2 

3 

101 

102 

103 

115 

123 

122 

121 

114 

55 

3 

3 

3 

4 

5 

16 

25 

24 

23 

15 

5 

4 

3 

5 

6 

7 

17 

27 

26 

25 

16 

5 

5 

3 

7 

8 

9 

18 

29 

28 

27 

17 

5 

6 

3 

9 

10 

11 

19 

31 

30 

29 

18 

5 

7 

3 

11 

12 

13 

20 

33 

32 

31 

19 

5 

8 

3 

21 

22 

23 

60 

103 

102 

101 

59 

50 

9 

3 

103 

104 

105 

116 

125 

124 

123 

115 

5 

10 

3 

105 

106 

107 

117 

127 

126 

125 

116 

5 

11 

3 

107 

108 

109 

118 

129 

128 

127 

117 

5 

12 

3 

109 

110 

111 

119 

131 

130 

129 

118 

5 

13 

3 

111 

112 

113 

120 

133 

132 

131 

119 

5 

14 

3 

23 

24 

25 

35 

46 

45 

44 

34 

1 

15 

3 

25 

63 

105 

92 

78 

62 

46 

35 

1 

16 

3 

105 

104 

103 

91 

76 

77 

78 

92 

1 

17 

3 

103 

60 

23 

34 

44 

61 

76 

91 

1 

18 

3 

25 

26 

27 

37 

49 

48 

47 

36 

1 

19 

3 

27 

66 

107 

94 

81 

65 

49 

37 

1 

20 

3 

107 

106 

105 

93 

79 

80 

81 

94 

1 

21 

3 

105 

63 

25 

36 

47 

64 

79 

93 

1 

22 

3 

27 

28 

29 

39 

52 

51 

50 

38 

1 

23 

3 

29 

69 

109 

96 

84 

68 

52 

39 

1 

24 

3 

109 

108 

107 

95 

82 

83 

84 

96 

1 

25 

3 

107 

66 

27 

38 

50 

67 

82 

95 

1 

26 

3 

29 

30 

31 

41 

55 

54 

53 

40 

1 

27 

3 

31 

72 

111 

98 

87 

71 

55 

41 

1 

28 

3 

111 

110 

109 

97 

85 

86 

87 

98 

1 

29 

3 

109 

69 

29 

40 

53 

70 

85 

97 

1 

30 

3 

31 

32 

33 

43 

58 

57 

56 

42 

1 

31 

3 

33 

75 

113 

100 

90 

74 

58 

43 

1 

32 

3 

113 

112 

111 

99 

88 

89 

90 

100 

1 

33 

3 

111 

72 

31 

42 

56 

73 

88 

99 

1 

1  -5.375  -1.375  1 

2  -5.1875  -1.375  1 

3  -5.0  -1.375  1 

4  -4.0  -1.375  1 

5  -3.0  -1.375  1 

6  -2.0  -1.375  1 

7  -1.0  -1.375  1 
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8  0.0 

-1.375 

1 

9  1.0 

-1.375 

1 

10  2.0 

-1.375 

1 

11  3.0 

-1.375 

1 

12  4.0 

-1.375 

1 

13  5.0 

-1.375 

1 

14  -5.375 

.  -1.1875  1 

15  -5.0 

-1.1875 

0 

16  -3.0 

-1.1875 

0 

17  -1.0 

-1.1875 

0 

18  1.0 

-1.1875 

0 

19  3.0 

-1.1875 

0 

20  5.0 

-1.1875 

0 

21  -5.375 

.  -1.0 

1 

22  -5.1875  -1.0 

0 

23  -5.0 

-1.0  0 

24  -4.0 

-1.0  0 

25  -3.0 

-1.0  0 

26  -2.0 

-1.0  0 

27  -1.0 

-1.0  0 

28  0.0 

-1.0  0 

29  1.0 

-1.0  0 

30  2.0 

-1.0  0 

31  3.0 

-1.0  0 

32  4.0 

O 

o 

1 

33  5.0 

o 

o 

1 

34  -4.65 

-0.65 

0 

35  -3.35 

-0.65 

0 

36  -2.65 

-0.65 

0 

37  -1.35 

-0.65 

0 

38  -0.65 

-0.65 

0 

39  0.65 

-0.65 

0 

40  1.35 

-0.65 

0 

41  2.65 

-0.65 

0 

42  3.35 

-0.65 

0 

43  4.65 

-0.65 

0 

44  -4.35355339059327373  -0.35355339059327373  1 

45  -4.0  -0.5  1 

46  -3.64644660940672627  -0.35355339059327373  1 

47  -2.35355339059327373  -0.35355339059327373  1 

48  -2.0  -0.5  1 

49  -1.64644660940672627  -0.35355339059327373  1 

50  -0.35355339059327373  -0.35355339059327373  1 

51  0.0  -0.5  1 

52  0.35355339059327373  -0.35355339059327373  1 

53  1.64644660940672627  -0.35355339059327373  1 

54  2.0  -0.5  1 
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55  2.35355339059327373  -0.35355339059327373  1 

56  3.64644660940672627  -0.35355339059327373  1 

57  4.0  -0.5  1 

58  4.35355339059327373  -0.35355339059327373  1 

59  -5.375  0.0  1 

60  -5.0  0.0  0 


61  -4.5 

0.0 

1 

62  -3.5 

0.0 

1 

63  -3.0 

0.0 

0 

64  -2.5 

0.0 

1 

65  -1.5 

0.0 

1 

66  -1.0 

0.0 

0 

67  -0.5 

0.0 

1 

68  0.5 

0.0 

1 

69  1.0 

0.0 

0 

70  1.5 

0.0 

1 

71  2.5 

0.0 

1 

72  3.0 

0.0 

0 

73  3.5 

0.0 

1 

74  4.5 

0.0 

1 

75  5.0 

0.0 

0 

76  -4.35355339059327373 

77  -4.0  0.5  1 

78  -3.64644660940672627 

79  -2.35355339059327373 

80  -2.0  0.5  1 

81  -1.64644660940672627 

82  -0.35355339059327373 

83  0.0  0.5  1 

84  0.35355339059327373 

85  1 . 64644660940672627 

86  2.0  0.5  1 

87  2 . 35355339059327373 

88  3.64644660940672627 

89  4.0  0.5  1 

90  4.35355339059327373 


91 

-4.65 

0.65 

0 

92 

-3.35 

0.65 

0 

93 

-2.65 

0,65 

0 

94 

-1.35 

0.65 

0 

95 

-0.65 

0.65 

0 

96 

0.65 

0.65 

0 

97 

1.35 

0.65 

0 

98 

2.65 

0.65 

0 

99 

3.35 

0.65 

0 

100 

4.65 

0.65 

( 

101  -5.375  1.0  1 


0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 

0.35355339059327373  1 
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102  -5.1875  1.0  0 


103  -5.0 

1.0  0 

104  -4.0 

1.0  0 

105  -3.0 

1.0  0 

106  -2.0 

1.0  0 

107  -1.0 

1.0  0 

108  0.0 

1.0  0 

109  1.0 

1.0  0 

110  2.0 

1.0  0 

111  3.0 

1.0  0 

112  4.0 

1.0  0 

113  5.0 

1.0  0 

114  -5.375 

.  1 . 1875  1 

115  -5.0 

1 . 1875 

0 

116  -3.0 

1 . 1875 

0 

117  -1.0 

1 . 1875 

0 

118  1.0 

1 . 1875 

0 

119  3.0 

1 . 1875 

0 

120  5.0 

1 . 1875 

0 

121  -5.375 

.  1.375 

►  1 

122  -5.1875  1.375  1 

123  -5.0 

1.375 

1 

124  -4.0 

1.375 

1 

125  -3.0 

1.375 

1 

126  -2.0 

1.375 

1 

127  -1.0 

1.375 

1 

128  0.0 

1.375 

1 

129  1.0 

1.375 

1 

130  2.0 

1.375 

1 

131  3.0 

1.375 

1 

132  4.0 

1.375 

1 

133  5.0 

1.375 

1 

13  3 
111 
111 

2  3  3 
1  1  1 
111 

3  16  3 

1111111111111111 
1  1  1 

4  16  3 

1111111111111111 

111 

5  16  3 

1111111111111111 

111 
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6  16  3 

1111111111111111 

111 

7  16  3 

1111111111111111 

111 

8  3  16 
111 

1111111111111111 

9  16  3 

1111111111111111 

111 

10  16  3 

1111111111111111 

111 

11  16  3 

1111111111111111 

111 

12  16  3 

1111111111111111 
1  1  1 

13  16  3 

1111111111111111 

111 

14  16  4 

1111111111111111 

1111 

15  16  4 

1111111111111111 

1111 

16  16  4 

1111111111111111 

1111 

17  16  4 

1111111111111111 

1111 

18  16  4 

1111111111111111 

1111 

19  16  4 

1111111111111111 

1111 

20  16  4 

1111111111111111 

1111 

21  16  4 

1111111111111111 
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L 


1111 

22  16  4 

1111111111111111 

1111 

23  16  4 

1111111111111111 

1111 

24  16  4 

1111111111111111 

1111 

25  16  4 

1111111111111111 

1111 

26  16  4 

1111111111111111 

1111 

27  16  4 

1111111111111111 

1111 

28  16  4 

1111111111111111 

1111 

29  16  4 

1111111111111111 

1111 

30  16  4 

1111111111111111 

1111 

31  16  4 

1111111111111111 

1111 

32  16  4 

1111111111111111 

1111 

33  16  4 

1111111111111111 

1111 
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B.3  Matlab  Plotting  Script 


function  [matrix.f ile]  =  gnu2mat(gnuf ile,  rows,  cols) 
y.  This  function  converts  a  GNU  file  generated  by  the  output  of  the 
y,  code  found  in  the  book  "Quick  Finite  Elements  for  Electromagnetic 
y,  Waves"  by  Pelosi,  Coccioli,  Selleri  to  a  matrix  format  readable 
y,  by  Matlab. 

y. 

y,  Syntax: 

y. 

y,  filename  =  gnu2mat(GNU_f ile,  #row,  #cols) ; 

y.  Perry  N.  Villanueva 
y.  1  December  1998 

y. 

index  =  1; 
for  i  =  Irrows 
for  j  =  l:cols 

matrix.f ile(i, j)  =  gnuf ile(index) ; 
index  =  index  +  1; 

end 

end 

figure 

setCgcf,  ’Position’,  [152  289  946  542]) 

subplot  211 

setCgca,  ’Units’,  ’inches’) 

'/.FOR  FIVE  CYLINDER  CONFIGURATION  (PBF)  ,  UNCOMMENT  NEXT  LINE 
set(gca,  ’Position’,  [0.75  3.5  9.55  2.0]) 

'/.FOR  ELEVEN  CYLINDER  CONFIGURATION,  UNCOMMENT  NEXT  LINE 
'/,  setCgca,  ’Position’,  [0.75  3.25  9.55  2.0]) 

contour (matrix_f ile,  8) 
grid 

axis  equal 

axis (  [0  cols  0  rows] ) 

'/,t  it  1  e  ( ’  Cont  our  Map  ’ ) 
xlabel ( ’ Pixels ’ ) 
ylabel ( ’ Pixels ’ ) 

subplot  212 

waterf all (matrix.f ile) 
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*/,colormap  (gray) 
eixis([0  cols  0  rows  0  3]) 

•/.setCgcf,  ’Position’,  [10  396  945  340]) 

set(gca,  ’Units’,  ’inches’) 

setCgca,  ’Position’,  [0.75  0.55  9.55  2.0]) 

setCgca,  ’Position’,  [0.8500  0.7500  9.5500  2.25000]) 

view(-4,72) 

*/,title( ’Field  Map’) 
xlabel ( ’ Pixels ’ ) 
ylabel ( ’ Pixels ’ ) 
rotate3d 

'/.figure 

'/,waterfall(matrix_f  ile) 


'/.axis([0  728  0  88  0  4]) 

'/.set(gcf,  ’Position’,  [10  396  945 

340]) 

'/,set(gca,  ’Units’,  ’inches’) 

'/,set(gca,  ’Position’,  [0.8500  0.7500 

9.5500 

2.25000]) 

'/.view  (-4, 65) 

'/.title (’Waterfall  Plot  of  Fields’) 

'/.xlabel  ( ’  Pixels  ’ ) 
'/.ylabel  (’Pixels’) 
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